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ABSTRACT 



We present a study of numerical behavior of a thickened flame used in Flame Cap- 
turing (FC, Khokhlov (1995)) for tracking thin physical flames in deflagration simu- 
lations. This technique, used extensively in astrophysics, utilizes artificial flame vari- 
able to evolve flame region, width of which is resolved in simulations, with physically 
motivated propagation speed. We develop a steady-state procedure for calibrating 
flame model used in FC, and test it against analytical results. Original flame model 
is properly calibrated with taking matter expansion into consideration and keeping 
artificial fiame width at predetermined value regardless of expansion. 

We observe numerical noises generated by original realization of the technique. 
Alternative artificial burning rates are discussed, which produce acceptably quiet 
flames (relative dispersion in propagation speed within 0.1% at physically interesting 
ratios of fuel and ash densities). 

Two new quiet models are calibrated to yield required "flame" speed and width, 
and further studied in 2D and 3D setting. Landau-Darricus type instabilities of 
the flames are observed. One model also shows significantly anisotropic propagation 
speed on the grid, both effects increasingly pronounced at larger matter expansion 
as a result of burning; these 2D/3D effects make that model unacceptable for use 
in type la supernova simulations at fuel densities below about 100 tons per cubic 
centimeter. Another model, first presented here, looks promising for use in flame 
capturing at fuel to ash density ratio of order 3 and below, the interval of most 
interest for astrophysical applications. No model was found to signiflcantly inhibit 
LD instability development at larger expansions without increasing flame width. The 
model we propose, "Model B", yields flames completely localized within a region 6 
cells wide at any expansion. 

We study Markstein effect (speed of the flame dependence on its curvature) for 
flame models described, through direct numerical simulations and through quasi- 
steady technique developed. By comparing results obtained by the 2 approaches we 
demonstrate that Markstein effect dominates instability effects on curved flame speed 
for Model B in 2D simulations for fuel to ash density ratio of about 2.5 and below. We 
flnd critical wavelength of LD instability by direct simulations of perturbed nearly 
planar flames; these agree with analytical predictions with Markstein number values 
found in this work. We conclude that the behavior of model B is well understood, 
and optimal for FC applications among all flame models proposed to date. 
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CHAPTER 1 
INTRODUCTION 

In the thesis we study features of propagating flame solutions in reaction-diffusion 
(RD) systems. Such flames arc ubiquitously observed in chemical combustion in 
terrestrial experiments and in industry, as well as in certain astrophysical phenom- 
ena, with fast composition transformation taking place in nuclear deflagration fronts 
propagating through stellar material. For most of the study we concentrate on a 
few specific RD models used in, or proposed here for use in simulations of deflagra- 
tion phase in type la supernova (SN la; plural SNe la) phenomenon. An important 
applied goal of the study is to propose an optimal flame model to be used in the 
simulations, the one that reliably reproduces expected physical propagation of a de- 
flagration front in SN la problem, without introducing unphysical numerical artifacts, 
such as numerical noises, excessive front instabilities and anisotropic behavior related 
to computation grid. 

We start with a brief overview of SN la phenomenon, providing motivation for 
the study, and discussing general physics of flames in this example system. We then 
proceed to describing techniques currently used for numerical simulations of SN la. 
The chapter is concluded with the list of problems the current simulations face that 
we address in this work, and with more detail on organization of the thesis. 

1.1 Supernovae: observations 

Extreme cases of variable stars, which are now classifled as novae and supernovae, 
attracted human attention since ancient times. The first accurate record of such 
an event, with documented dates (185 AD) and position (close to a Centauri) is 
due to Chinese astronomers (Zhao et al. (2006)). Crab Nebula Ml and the famous 
short-period pulsar within are remnants of another supernova, observed by Chinese, 
Japanese and Arab astronomers in 1054; for 23 days it was bright enough to be visible 
in daylight, and at nights it was visible with naked eye for 653 days. In 1572 Tycho 
Brahe measured parallax of another star of this class (to find no detectable paral- 
lax); he also recorded its dimming with time. That was the first well-documented 
supernova in western world; Brahe's book, "Be Stella Nova", gave rise to the name 
we currently use for a different explosive phenomenon, much more frequent and dim 
than supernovae. 

There were more than a hundred classical novae (now understood as a result of 
explosive burning of accreted hydrogen-helium fuel on a surface of a White Dwarf 
star in binary systems, when certain critical mass of the fuel is accumulated) studied 
by 1919, when Lundmark suggested that the distance to M31, "Andromeda Neb- 
ula", was 7x 10^ light years (at that time most astronomers believed that everything 
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observed on the sky belonged to our Galaxy). Even with this significantly underes- 
timated distance it followed that S Andromedae explosion of 1885 was extragalactic, 
and hence was more luminous than classical novae studied by about 3 orders of mag- 
nitude (Lundmark (1920)). With improved distance measurements Baade & Zwicky 
(1934) suggested to differentiate between novae and supernovae, in the way as used 
nowadays. 

Typical supernova brightness rises rapidly for a few weeks, reaches "maximum 
light" , then drops within months, approximately exponentially after a few weeks of 
less regular transient dimming. At maximum light SNe are among the brightest 
objects in the Universe, often shining brighter than the whole host galaxy. Zwicky 
(1938) proposed collapse of a star made of ordinary matter into a neutron star as 
a mechanism for SN phenomenon; currently accepted scenario of certain types of 
core-collapse supernovae is based on that idea. All supernovae are understood now 
to be a result of cataclysmic event happening to some stars in late stages of their 
evolution, after almost all hydrogen and helium are transformed into carbon and 
heavier elements. As observed energy emitted within a few weeks following such 
an event (around maximum light), is comparable to energy produced by the Sun 
over its multi-billion year history, it is clear that SN event must involve drastic 
transformation of the whole star structure. In fact, the transformation of most of 
the star composition takes just a few seconds based on current understanding of the 
physics of SNe. 

Historically supernovae are categorized into a few types based on their spectral 
features. Supernovae of type II show hydrogen absorption lines in the spectrum 
near maximum luminosity, those of type I show no hydrogen. Supernovae of type I 
showing strong Si II lines belong to class la; those with no silicon lines are further 
subdivided into types lb and Ic if they show or lack helium lines, respectively. 

Observationally SNe la are special in that vast majority of them show close peak 
luminosity and light curve decline features, whereas variety within other types is 
significant. Typical SN la light curve (brightness vs time) shows initial rise interval 
about 20 days long, followed by about as fast decline, about 3 magnitudes within a 
few weeks, then dimming more slowly by about a magnitude per month. For most 
spectroscopically normal SNe la the absolute bolometric magnitude at maximum 
light lie in narrow interval from —19.4™ to —18.7"^. Phillips (1993) observed that 
this scatter in peak luminosities is correlated with light curve width (characterized 
by Ami5(B), star magnitude drop within 15 days after maximum light); this Phillips 
relation makes it possible to accurately predict SN la absolute magnitude based on 
observed dimming rate. Coupled with extreme brightness of supernovae this makes 
SN la excellent standard candles for probing cosmological distances. Pecuharities in 
observed SN la magnitude as a function of their redshift was a first indication that 
the Universe is expanded with (positive) acceleration (see Leibundgut (2001) for a 
review) . 
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SNe la also differ from SNe of different types in their environment: la is the 
only supernova type observed in elliptical galaxies, although they still occur about 
twice more frequently in younger stellar populations in spiral and irregular galaxies 
(Cappellaro et al. (1997)). 

These features, together with low luminosities of SN la in radio. X-ray and gamma 
bands suggest less massive progenitor of SNe la than those of other types of super- 
novae, and no collapse happening in the process of explosion. All other supernova 
types are believed to involve core collapse to a neutron star or a black hole as a main 
source of energy powering disruption of and energy deposition into the envelope. 

1.2 SNe la: physical models 

Below we briefly overview physical models of SN la, mostly paying attention to 
physical issues relevant for our study. See, e.g., Hillebrandt & Niemeyer (2000) for 
more detail on explosion models. 

1.2.1 Overview 

SNe la are believed to be a result of thermonuclear explosion of a White Dwarf (Hoyle 
& Fowler (1960), Arnett (1969)), a late stage of evolution of hght stars, with main 
sequence mass below about 8 solar masses. All models involve heavy White Dwarfs 
(WDs), with masses close to or exceeding solar mass Mq. Such WDs are formed from 
contracted central parts of stars after hydrogen and helium burning there is complete 
and H/He-rich shell is dissipated in space. They are composed of -"^^C and -"^^O mostly, 

and are supported by the cold pressure of degenerate electrons. WDs with smaller 

on 0/1 

amounts of carbon, but rich in heavier elements, Ne and Mg, arc also candidates 
to undergo a thermonuclear explosion, however some calculations (Gutierrez et al. 
(1996)) indicate that core collapse is a probable outcome of near-central ignition in 
such WDs. 

For all specific models (briefly discussed below) not overtly contradicting observa- 
tions successful result of the explosion is stellar material transformed mostly into ^^Ni 
with considerable amount of intermediate-mass elements (Si, S, Ca) within about 3 s 
of proper explosion; enough energy is released for total disruption of the WD. Explo- 
sion itself is practically invisible due to opacity of dense hot ashes of the explosion; 
signiflcant part of the released energy is spent to overcome gravitational potential of 
the compact star, with excess transformed into kinetic energy of expanding material. 
A few days after explosion, expanding envelope starts looking bright to an external 
observer, though it remains optically thick until much later times. Double /3-decay of 
^^Ni to ^^Co {ti/2 = 6.08 d) and then to stable ^^Fe {ti/2 = 77.2 d) powers the light 
curve (Truran et al. (1967), Colgate & McKee (1969)), that is it heats the expanding 
gas, which captures emitted in the decay positrons and 7— quanta and eventually 
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reemits light outside, mostly in optical band. Both the initial light-curve decline rate 
(governed by slow diffusion of light out of the gas nebula) and peak luminosity (Ar- 
nett (1982)) turn out to be proportional to amount of ^^Ni produced, thus providing 
an explanation for Phillips relation. 

Specific models proposed differ in particulars of progenitor, ignition and burning 
mechanism. The most natural model to account for uniformity of most SN la events, 
single-degenerate scenario with nearly Chandrasekhar mass M(j]^ (the maximal mass 
of the WD that degenerate electron pressure can support, Chandrasekhar (1931); 
about 1.4 M0), involves ignition near the center of a heavy WD as a result of ther- 
monuclear runaway, readily happening in highly-degenerate fuel (C+0) in central 
part of a WD with near-critical mass. The progenitor WD is believed to be formed 
with lower mass, and to have its mass increased by accretion from a companion giant 
or main-sequence star. Parameters of the close binary system must be rather special 
so that accreting hydrogen or hehum steadily burned into C/0 on the surface of a 
WD, without nova outbursts (which do not on average lead to increase in a WD 
mass); cf. Nomoto et al. (2000) for study of the accretion process. 

The amount of nickel produced, and thus the maximal luminosity (if explosion 
does unbind the WD) is likely to depend on WD composition, rotation, and on 
specifics of explosion process, thus leading to variation in nickel mass produced. In 
the model of Arnett (1969) explosion starts as detonation, and it is in this detonation 
front sweeping through entire star where all "fuel" , C+0, is transformed into ash; the 
latter consists almost completely of nickel after detonation at original high density 
of WD material. This contradicts observations, which show a few tenths of solar 
mass of intermediate mass elements produced by all SNe la. Pure detonation model 
(in near- Chandrasekhar mass WD) is thus ruled out by observations (Arnett (1974), 
Ostriker et al. (1974)). It is believed that "explosion" in a single degenerate near- 
Mqyi model is started as defiagration, possibly turning into detonation at a later 
stage. We will discuss this in more detail in the next section. 

Other class of models involve WDs with masses substantially below Mq^^. While 
such lighter WDs (still massive enough to consist mostly of carbon/oxygen) will not 
self-ignite, detonation front may be sustained once initiated, according to simulations. 

Woosley & Weaver (1994), Livne & Arnett (1995) present ID and 2D simulations, 
in which detonation is initiated on the surface of C/0 subchandrasekhar WD, by 
self-detonating (through runaway at the bottom) accreted helium shell of ~ O.2M0. 
Detonations at lower densities in such lighter WDs can produce significant amounts of 
intermediate-mass elements, though some features observed in simulations (like large 
amounts of ^^Ti produced) are not supported by observations. Models of this type are 
attractive as binary systems containing lower mass WDs are more abundant; however 
accumulation of significant mass of He is still required for the mechanism to work; 
also, due to a range of possible progenitor masses there is no natural explanation for 
uniformity of SN la peak luminosities. Finally, simulations performed in 2D with 
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postulated axial symmetry are just not reliable, as 3D simulations of the last decade 
suggest. 

In double degenerate scenario (Webbink (1984), Iben & Tutukov (1984)) igni- 
tion is initiated by collision of 2 sub-MQi^ WDs. This scenario requires very close 
binary system of WDs orbiting around each other: their separation must not exceed 
~ 10^ km so that the merging happened on cosmological timescales, assuming grav- 
itational radiation as the main damping mechanism. Like the model in the previous 
paragraph there is no special mass for merging system, thus significant dispersion 
in maximum luminosities is to be expected. Such non-standard mechanisms may be 
an explanation for a small number of SNe la with peculiar spectra and substantially 
high or low peak luminosity. 

1.2.2 Deflagration and detonation in a White Dwarf: physics 

As reviewed in the previous section, nuclear explosion of a near-MQ^ WD is currently 
considered the best explanation for SN la. What is usually loosely referred to as "ex- 
plosion" physically corresponds to one of 2 possible modes of fast fuel transformation 
(through a network of thermonuclear reactions) into ash, in thin front propagating 
through the system (WD). Reaction rates show sharp dependence on temperature, 
of Arrhenius type. Temperature in the reaction zone (within the front) must be 
raised to approximately reaction activation energy for the reaction to progress with 
substantial rate. 

Detonation mode is characterized by supersonic front speed; adiabatic compres- 
sion by significant pressure jump across the front is the mechanism raising the tem- 
perature high enough for reactions to proceed. Heat released in the reactions sustains 
the shock wave propagation. Due to fast detonation front velocity, if detonation is 
initiated near WD center as the start of the explosion, the whole WD will be swept 
by the front at nearly original high density — density distribution has just no time 
to change. Reaction goes all the way through nuclear statistical equilibrium (NSE) 
composition at densities above about 5 x lO^g cm""^; and this composition is mostly 
iron-group elements (having the highest binding energies per nucleon) at high fuel 
densities. This is why too little lighter elements is produced in pure detonation 
model. 

In defiagration mode a fiame (thin zone with large gradients of temperature, 
density, composition; it contains reaction zone) propagates into the fuel by diffusion- 
reaction mechanism. Temperature is raised in front of the flame due to heat con- 
duction from the hot reaction zone. Temperature distribution has characteristic for 
diffusion processes exponential tail deep into cold fuel; right in front of the reaction 
zone (where still the fuel is not depleted in the reaction) the temperature continu- 
ously approaches high flame temperature, leading to high reaction rate; fiame this 
way propagates forward. Thermal energy released in reaction zone prevents temper- 
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ature profile from smearing out by lieat conduction; steady solutions (of the form 
of a kink wave) typically exist describing propagating flame with constant in time 
distribution of physical quantities within, in the flame rest-frame. See Chapter 2 for 
details; Fig. 2.1 shows ash fraction distribution in one particular flame model. 

Species diffusion, and their production/depletion in the reaction zone contribute 
to distribution of composition and temperature near the flame (see Sec. 5.1 below 
for details). In terrestrial combustion in nearly ideal gases burning is affected to 
approximately the same degree by diffusion as it is by heat conductivity, as Lewis 
number (ratio of heat diffusivity to matter composition diffusivity) Le 1. For WD 
matter thermal conductivity by far dominates other transfer phenomena, Le ~ 10^, 
thus species diffusion plays negligible role in flame propagation. 

Laminar flame speed as estimated in Timmes & Woosley (1992) depends on fuel 
composition (that is local WD composition in front of the flame) and density; it 
is about 80 kms"-*^ in central region of 0.5C-I-0.5O near-M^h WD (central density 
Pfuel ~ 2 X 10^ g cm""^) and lower at smaller densities in outer layers. This velocity 
can be roughly estimated by equating characteristic energy release timescale, 

ireac ~ 1/-R 

with heat diffusion timescale (assuming comparable or smaller contribution to flame 
dynamics from species diffusion, i.e. Le ~ 1 or greater) 

^diff ~ ■ 

K 

R here stands for some characteristic value of reaction rate in flame zone, such that 
the average energy release rate in the flame (dQ/dt) = Rq, q being total heat release 
of the reaction; k, = \/{cyp) is temperature diffusivity coefficient (A meaning heat 
conductivity, p density and cy specific heat). This balance of timescales determines 
flame width W , which further provides rough estimate for laminar flame speed: 

W « y/K/R, Si W/treac ~ VkR . 

See Zeldovich et al. (1985); WiUiams (1985) for more accurate estimates, and for 
more on physics of deflagration. 

1.2.3 Pure deflagration models 

Whether deflagration alone can account for observable features of SNe la is not 
straightforward to answer. Laminar flame speed is far less than sound speed any- 
where in the WD, thus flame propagation is nearly isobaric; star structure changes 
as heat is released in the flame, WD expands, this hydrodynamic evolution must be 
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computed simultaneously, coupled to flame propagation calculations. Early ID simu- 
lations suggested that fast WD expansion after deflagration sets in quenches burning 
too early, so that not enough energy is released to unbind the star; WD contracts 
back after burning stops. It was understood, however, that in reality the flame will 
never propagate as an ideal spherical flame with its laminar burning speed, due to 
instabilities discussed below. 

It is known that deflagration fronts are subject to a number of instabilities, which 
can drastically affect flame propagation. As soon as transverse flame dimensions ex- 
ceed characteristic instabihty scales Acr flame surface starts developing various fea- 
tures (bubbles, wrinkles — the geometry depends on speciflc instability involved) 
increasing its surface compared to a smooth surface stable laminar flame would 
have had; this effectively increases integral burning rate. In flamelet regime, that 
is when characteristic spatial scale of flame surface perturbations by flame instabil- 
ities and background turbulence signiflcantly exceeds laminar flame width, kinetics 
of the burning (and hence distribution of physical quantities) within the flame is not 
affected significantly; integral burning rate is increased approximately in proportion 
to flame surface increase (compared to stable laminar burning). When a flame is 
perturbed on scales comparable to or less than its laminar width burning description 
based on laminar flame physics becomes not applicable, the flame is generally torn 
apart; such regime is called distributed burning (Damkohler (1940)). 

For nuclear deflagration in degenerate WD matter, where kinematic viscosity 
and species diffusivity are negligible compared to temperature diffusivity k, hydro- 
dynamic instabihty of the type (abbreviated LD below) studied by Landau (1944), 
Darrieus (1938) is important. Besides, the flame is subject to Rayleigh-Taylor (RT) 
instability (basically a buoyancy instability — occurring when lower density reaction 
products support denser cold fuel above in WD gravitational field, directed towards 
its center), and to Kelvin-Helmholtz (KH) instability, perturbing shear fiow inter- 
faces, hke on sides of rising RT bubbles. LD instabihty is most important on scales 
comparable with flame width. Such scales are never resolved in full star simulations 
(flame width is submillimeter for most of deflagration phase Timmes & Woosley 
(1992)). On the largest scales RT instability dominates, driving hydrodynamics on 
the scale of the star. See Khokhlov (1994, 1995) for features of RT-driven burning in 
cylinders with uniform gravitational fleld, and in 3D SN la simulations. One result 
of studies in cylindrical geometry important in applications (see next section) is that 
after certain transient time steady (in average) burning sets in uniform gravity driven 
RT burning, with integral burning velocity ("turbulent flame speed", St) depending 
only on geometry of the cylinder cross-section, and independent of laminar flame 

speed. 

St ~ (1.1) 

is the value conflrmed in Khokhlov (1995), Zhang et al. (2007) as reasonably accurate. 
A here stands for linear scale of cylinder cross-section (in the cited works horizontal 
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cross-section was a square with side A); A denotes Atwood number for fuel-ash, 

^ = (Pfuel - Pash)/ (Pfuel + Pash), (1-2) 

g stands for gravitational acceleration. 

The smaller the laminar speed the more intricate surface of the flame develops, its 
area in steady regime scaling inversely proportional to Si, laminar speed of the flame. 
If Si is increased, on the other hand, the surface becomes smoother; when Si > Sf 
(as given in (1.1)) flame is essentially flat and propagates with its laminar speed; 
that is for any geometry of the system flames with large enough Si are stable to RT. 
See Khokhlov (1995) for more on self- regulation mechanism behind this behavior; 
qualitatively, for flames with larger laminar speed laminar burning overruns slower 
developing smaller RT bubbles, thus only larger-scale bubbles (which rise faster) 
develop — the surface is polished by burning on lower scales. 

In non-steady 3D setting in deflagrating WD characteristic large scale burning 
speed St on large scales increases as the flame gets farther from the WD center, as 
Atwood number, gravity and available tangential dimension increase; in addition, the 
instability had been allowed more time to develop before the flame got to such larger 
radius. Laminar flame speed, on the other hand, rapidly decreases, as fuel density 
decreases. Hence, early ID simulations grossly underestimated total burning rate, 
except at the very center of the WD. ID models are used for certain purposes till 
now, however the flame speed is prescribed based on estimates of Sf from 3D analysis. 
Often this burning rate is being left as a (time dependent) free parameter, tuned to 
get desired features on the output of the simulation, close to observations; some 
calculations of this type are notably successful in yielding results closely resembling 
observed SN la features (Nomoto et al. (1984)). 

Real 3D simulations are needed to judge about viability of a pure deflagration 
model. These still use artificial techniques to propagate the flame with right velocity 
(improving one of such techniques, proposed in Khokhlov (1994), is a subject of this 
thesis) as the scales physically governing flame propagation are not resolved in any 
WD-scale 3D simulations. However the techniques aim to get the flame propagating 
with correct, physically motivated effective speed to reproduce, on resolved scales, 
speed of real thin physical flame; no tuning of flame speed with the purpose of getting 
particular output results is used. We will discuss the techniques used in the next 
section; more detail on implementation and the results may be found in Khokhlov 
(2000), Gamezo et al. (2003), Gamezo et al. (2005), Reinecke et al. (2002a), Niemeyer 
et al. (2002). 

No consensus exists right now whether deflagration alone may explain a normal 

type la supernova. With currently used subgrid models, prescribing flame speed so 
that to get total burning rate as the physical one (as the latter is currently under- 
stood), of submillimeter-thick nuclear flame wrinkled on unresolved scales, different 
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groups get enough energy to unbind the WD, for near-central single-spot ignition. 
The amount of energy released in such deflagration models is still lower than in 
normal SNe la, as is amount of iron-group elements produced. Significant amount 
of unburned carbon and oxygen remains close to the center of the star, a result of 
sinking RT fingers of cold fuel. After homologous expansion of unbound WD sets in 
(about 10 seconds after original ignition) these fingers arc observed as large amounts 
of unburned material expanding with low velocities; on the other hand in real SNe la 
light elements are observed almost exclusively in outer layers of expanding nebula, 
that is having largest expansion speeds. 

Attempts were undertaken (Ropke & Hillebrandt (2005)) to fix some of these 
shortcomings with extending burning simulation at lower densities, in distributed 
burning regime (starting at densities below about 10^ g cm""^, when flame preheating 
zone becomes wider than Gibson scale. 

Another (effective) approach (Travaglio et al. (2004); Ropke et al. (2007)) is 
based on multipoint ignition. Igniting at several points distributed in the central 
part of a WD simultaneously reduces the amount of unburnt fuel near the center 
and increases the amount of nickel produced. Say, simultaneous ignition in 1600 
spherical spots in Ropke et al. (2007) resulted in asymptotic kinetic energy of ejecta 
of 8.1 X 10^^ erg, with 0.606 of iron group elements produced, including 0.33 
of ^^Ni. Corresponding results in Reinecke et al. (2002b) with single spot central 
ignition (albeit smaller resolution) are 4.5 X 10^0 erg, 0.5 M© and 0.3 M©. Simulation 
in Travaglio et al. (2004), in one octant of a full WD (octahedral symmetry assumed) 
with 30 bubbles (larger than in Ropke et al. (2007), 10 km radius vs 2.6 km; and 
located closer to WD center) ignited at t = resulted in 6.5 x 10^^ erg asymptotic 
kinetic energy of ejecta and 0.418 Mq of ^^Ni produced. Kinetic energy of ejecta of 
a typical observed SN la is slightly above 10"*^ erg. Exact pre-explosion conditions 
(smoldering phase) at the center of a WD are not known well, so it is not obvious if 
single point ignition is the best choice of initial conditions. Yet significant dependence 
on this choice of initial conditions is noteworthy, even though dispersion in Ni masses 
quoted above would lead to dispersion in brightness within 0.8"^, comparable to 
dispersion among observed la's. 

1.2.4 Delayed detonation models 

The best fit with observations is obtained in hybrid models, where "explosion" starts 
as deflagration, and later on detonation is triggered in some spot of the WD pre- 



1. This is the scale on which average turbulent velocity fluctuations are equal to the laminar flame 
speed. On lower scales the turbulent fluctuations are smaller, therefore flame propagation is insen- 
sitive to turbulence on such lower scales. Above Gibson scale flame surface is bent and stretched by 
background turbulent flows. The turbulence on all scales is mostly generated (through Kolmogorov 
cascade) by large-scale RT motions, augmented with vorticity generated by KH instability. 
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expanded during the deflagration phase. Detonation happening at lower densities 
under such a scenario produces sufficient amount of intermediate mass elements (in 
outer layers of preexpanded WD), transforms RT fingers of fuel near WD center into 
iron-group elements (together with the whole dense central part: detonation front 
is smooth, on scales characteristic for RT instability, in contrast with intricately 
corrugated deflagration front surface), leaves less fuel unburnt and produces more 
energy. Exact mechanism, what triggers detonation, is not yet agreed upon. When 
originally proposed, Khokhlov (1991) assumed deflagration to detonation transition 
(DDT) at smaller densities (a few times 10^ gcm~"^, characteristic for distributed 
burning regime) as the most probable mechanism. DDT is observed in terrestrial 
chemical flames, however it is more problematic to explain it in unconstrained WD 
setting. 

Different indirect mechanisms for triggering detonation were proposed over time, 
involving creation of extended regions with large temperature gradients by colliding 
masses of WD material with different velocities. Pulsating detonation scenario in- 
volves recollapse of the WD after deflagration phase in which not enough energy was 
released to unbind the star. The original implementation (see Hoeflich & Khokhlov 
(1996)) seems not plausible now, as accurately modeled 3D deflagration with central 
ignition does ordinarily produce enough energy to unbind a C-l-0 WD, as reviewed 
in the previous subsection. 

Asymmetric variations of this scenario, with slightly offcenter ignition, Gravitati- 
onally-Conflned Detonation (GCD, Plewa et al. (2004)), Detonating Failed Deflagra- 
tion Model (DFD, Plewa (2007)) do trigger detonations in 3D simulations. In these 
only a small fraction of fuel is burnt before the bubble driven by buoyancy breaks 
through WD surface, creating surface wave intensive enough to trigger detonation 
after colliding at diametrically opposite point of the WD surface. The main problem 
of these models is to get enough preexpansion before detonation so that signiflcant 
amount of intermediate mass elements be produced, to agree with observations. Pre- 
expansion, depending on the mass burnt before the bubble breakthrough, strongly 
depends on the position of original ignition site. The problem is thus to propose 
some robust mechanism to rule out signiflcantly offcenter ignitions. These models 
tend to overproduce iron-group elements in all simulations to date. 

1.3 Simulations of deflagration phase 

Here, after a brief overview of standard hydrodynamics, we describe ff am c- capturing 
(PC) technique (Khokhlov (1995)) used in SN la simulations for tracking ffame po- 
sition and prescribing heat release in artificially thickened (to be resolved in sim- 
ulations) fiame zone with intent to reproduce, on resolved scales, results of real, 
non-discretized, hydrodynamics. We then describe some problems of the original 
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implementation of FC found in this study, and outline topics presented in the thesis, 
mostly related to improving the FC model. 

1.3.1 Hydrodynamics involved in simulations 

Deflagration is described by a standard set of hydrodynamic equations with reaction 
source terms. These include a mass continuity equation; Navier-Stokes equation for 
velocity of the matter subject to gravity, pressure gradients and viscosity; transport 
equation for internal energy, with a source term describing heat release in reactions 
taking place in the system; and finally equations describing diffusion and reactions 
for species involved in reactions. Equations of state, expressing pressure and internal 
energy as functions of density and temperature make the system of equations closed. 
For SN la problem the system is usually simplified by neglecting viscosity and dif- 
fusion of species (as thermal conductivity by several orders of magnitude dominates 
other transfer phenomena: Prandtl number is small, Lewis number is large). 

For large scale simulations the system is modified further: detailed reaction net- 
works are not used, instead only a few species are taken into account, with model 
reaction rates for these, found separately in such a way that the reduced system 
imitated kinetics of the full system with acceptable accuracy. Such a simplification 
is used in chemical combustion simulations as well: this makes it possible to save 
computational resources, to be able to use higher resolution for example, when it 
is more important to resolve some critical hydrodynamic scales then to get exact 
distributions of species. This is usually the case in engineering applications, when, 
for example, knowing pressure or temperature distribution in combustion chamber 
is more important than knowing detailed chemistry of burning process^. 

In SN la problem using detailed reaction network in large-scale hydrodynamic 
simulations makes little sense in principle, as physical reaction zone (across which 
at least some species concentrations change significantly) is never resolved in such 
simulations, and this situation will not change in any foreseeable future. As the total 
simulation box size must be of the WD size, that is a few thousand kilometers, grid 
cells are of order 1 km size for most detailed of current simulations. Flame width is 
submillimeter for most of the deflagration (at density in front of the flame Pf^gl ~ 
lO^g cm""^ and above, Timmes & Woosley (1992)). Because of this disparity, all full 
star simulations evolve flame surface indirectly, without resolving flame width"^. 



2. The latter is studied separately, for calibrating the reduced scheme in part, in simplified 
(usually steady ID) setting, allowing to concentrate on kinetics rather than (trivial in such setup) 
aerodynamics. 

3. The flame propagates directly on its own when the sot of hydrodynamic equations described 
above is solved exactly. Then, as described in Sec. 1.2.2, interplay between heat release in the flame 
and heat conduction into cold material in front of the flame leads to flame propagating forward, 
with correct (by definition) speed. In simulations described in this paragraph, however, derivatives 
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One approach {level set method, LSM; used by Munich group in part, Reinecke 
et al. (2002a); Ropke et ah (2007)), motivated by tiny width of real flames, is to 
(formally) describe the flame as infinitely thin surface, on which matter density, ve- 
locity and other physical variables have hydrodynamically prescribed discontinuities. 
One then propagates this discontinuity surface with the speed prescribed so that 
to mimic, in average, real flame region evolution. This surface is represented as 
the zero-level of certain "level function" G: flame manifold|^ = {r{t) '■ G{r,t) = 0}. 
Time evolution of G (advection plus normal propagation into the fuel with prescribed 
"flame" speed) used strives to ensure the needed evolution of the flame surface. Var- 
ious numerical tricks are employed to take care of the tendency of G to develop 
unphysical peculiarities. Hydrodynamical solvers (based on piecewise-parabolic re- 
construction. Woodward & Colella (1984)) used in conjunction with this formally 
inflnitely thin flame model spread discontinuities over several grid spacings, exact 
fractions of fuel/ash in grid cells intersected by the "flame surface" have more of 
numerical significance for the scheme than physical meaning. 

We will stick to another fiame-tracking prescription (described below), where 
the deflagration front is represented with artiflcial thick "flame" . Quotation marks 
are used below to distinguish the artiflcial numerical flame wherever confusion with 
real physical flame seems possible. We also refer to this numerical construction 
as "flame zone" as it is intended to be evolved in such a way so that to contain 
the physical flame within. The "flame" is governed by the same reaction-diffusion 
equations as physical flames, which makes it clear what behavior of the "flame" to 
expect (instabilities, interaction with background turbulence and sound waves), no 
new physics is introduced. This seems a clear advantage over LSM. 

1.3.2 Flame capturing 

Flame Capturing (FC) technique (Khokhlov (1995)) employs artiflcial scalar quantity 
/(r, t) e [0; 1], "reaction progress variable", to track flame evolution: / = in the 
unburnt material, and changes to 1 behind the "flame" , when the burning is complete. 
The "flame" in the sense of this numerical technique is a region where / is neither 
nor 1 (but strictly between); it is made a few grid sizes thick by appropriately choosing 
parameters governing / evolution (see below^); the value of 1 — / is intended to mimic 
real physical fuel fraction in fuel/product mixture within the grid cell in the flame 
region. 



are modeled via finite differences of quantities in adjacent grid cells. This modeling is designed to 
be accurate enough for fluid motions on scales exceeding a few grid spacings, yet obviously such 
differences cannot model derivatives of quantities changing significantly on subgrid scales, which is 
the case for reaction species abundances. 

4. This in essence is a particular realization of artificial viscosity approach, Neumann & Richt- 
myer (1950). 
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/ is evolved via a diffusion-type equation, 

^ + u-V/ = V(i^V/) + $(/). (1.3) 

Pfiysical quantities (pressure, temperature, composition, matter velocity) are simul- 
taneously evolved through a standard system of Euler equation (5.2), diffusion-type 
equations with reaction terms for internal energy density (heat conduction equation) 
and species concentrations (if tracked separately; in standard FC realization they 
are just set to be a linear function of /) , and equations of state (in code ALLA that 
we use for most non-steady simulations these are implemented as functions taking 
matter density and internal energy density as input, and returning pressure and 
temperature). / is coupled with physical variables through advection (by local mat- 
ter velocity u in (1.3)); and through heat-release term in hydrodynamic equations, 
which is governed directly by /: 6Q = qdf, that is heat is released hnearly with / 
increasing (reaction progressing) up to q, heat release of complete nuclear burning, 
per unit mass. (This q depends on local fuel composition and density.) Artificial 
diffusivity K and artificial "burning rate" $(/) were prescribed (Khokhlov (1995)^) 
so that to make the "flame" ~ 4 A thick, A denoting grid zone size, and propagating 
with prescribed velocity. That original prescription in effect led to systematic error 
in flame velocity, due to matter expansion (in the process of burning) neglected when 
estimating front speed in a system (1.3). This was corrected in Zhiglo (2007) (and 
will be described in the next chapter). 

Scales larger than the flame width, which govern flame instabilities development 
(LD and RT), are not resolved in SN la simulations as well — for most of the simu- 
lation, at large densities; a portion of a flame within any grid cell is thus not smooth, 
but corrugated by instabilities. Real flame surface therefore exceeds the surface im- 
plied from numerical observations on grid scale. This is an extra complication, not 
encountered in shock-tracking techniques used similarly to FC in certain simulations 
with shocks. To address this issue all simulations include some subgrid model for 
prescribing renormalized, corrected for missing scales velocity, larger than physical 
laminar flame speed. 

Observations of stationary RT-accelerated burning in Khokhlov (1995) led to 
using "turbulent flame speed" St, defined through (1.1) where for A the grid spacing 
is used, as the value for artificial "flame" propagation speed ensuring correct integral 
burning rate. When the laminar flame speed Si exceeds Sf, before developed RT- 
driven burning regime sets in, the "flame" is required to propagate with laminar speed 
Si- See Zhang et al. (2007) for further verification of this flame speed prescription, 
and for proposed technique modification in transient, non-stationary setting. We will 



5. Only constant K was used before Zhiglo (2007). Eq. 1.3 here is written in the form we use 
below with more general flame models, in which is a function of /. 
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not touch below on the best prescription for the "flame" speed, but concentrate on 
the very model propagating the "flame" with a given speed. 

1.3.3 Goals and scope of the thesis 

To hope to get meaningful physical results in simulations it is critically important 
to understand the properties of the model used for tracking the flame region and 
for estimating heat release there. It is important that flame tracking model itself, 
in part, did not introduce unphysical artifacts in uncontrolled way into simulations. 
One needs to study model behavior in density range between ~ 2 x 10^ g cm~*^ (cen- 
tral WD density) and a few xlO^ g cm~*^ (when distributed burning sets in, and 
hypothetical deflagration to detonation transition could take place), and favor mod- 
els not introducing much noise into simulations, and not demonstrating significant 
instabilities of their own, unless there is a hope that the latter mimic instabilities of 
actual physical flame zone. 

When this work was started no analysis of features of the flame model used 
for FC existed, apart from initial study of isotropy of "flame" propagation on the 
grid in Khokhlov (1995), together with analytical estimates that led to prescription 
for parameters in (1.3); the parameter values proposed there were used essentially 
verbatim by several groups thereafter (Gamezo et al. (2003); Plewa et al. (2004); 
Brown et al. (2005)). By now some results on numerical noises generated by flame 
models, including description of a new burning rate proposed for use in FC, were 
presented by our colleagues (Jordan et al. (2008)). 

Numerical simulations always pose a question how reliable the results are, what 
in the results is due to actual physics of the system studied, and what is an artifact 
of the numerical scheme/approximation/model used. This question is particularly 
important for systems so complex that simulations are the only way to get an es- 
timate of the results. With these concerns results are always tested against other 
simulations, with different resolution, utilizing different discretization algorithm or 
based on analytically different solver. A suite of problems with known analytical 
behavior is used to compare these known results with results of simulations, albeit 
on simple systems. With all this care every so often discrepancies between results of 
different groups appear. Khokhlov (1995) reports that 3-cell wide cylindrical flame 
propagates isotropically on a square grid. Same FC model implemented by Niemeyer 
required 8-cell wide flame to make anisotropies acceptable. CCD model reach robust 
detonation conditions in 3D in FLASH group simulations (Jordan et al. (2008)), the 
same model never detonates in simulations by Ropke et al. (2007). SN la simulations 
in octant with central ignition show clear tendency of the flame to preferentially de- 
velop features along grid directions from the very beginning of the simulation, which 
is a clear numerical artifact (this further facilitates fast RT instability development). 
Concerns like this were a part of motivation to study flame model features in detail. 
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In Chap. 2 we present results obtained mostly without using actual non-steady 
hydrodynamical solvers. We describe a method we developed for finding fiame speed 
by solving eigenvalue steady-state problem. We use that method to calibrate the 
model flame used in Khokhlov (1995) for flame capturing. In original calibration 
in Khokhlov (1995) matter expansion (we characterize it with expansion parameter 
^ — Pash/ (Pfuel ~ Pash) ) ^ result of burning was neglected; we correct for this and 
propose additional improvement, allowing to keep fiame width independent of matter 
expansion (this expansion increases with fuel density decreasing, due to decreasing 
degeneracy of the electron gas) . We present analytical solution for the original model 
fiame profiles, which allows us to directly check the eigenvalue numerical scheme. 
Flame profile (distribution of / in space for steadily propagating fiame solution of 
(1.3)) has an exponential tail; we propose another class of models, with /-dependent 
diffusivity coefficients, that produce flames localized in space (having finite total 
width, with no tails). We argue about the advantages of models having such finite 
flames, and find certain representative model of that class with nearly linear fiame 
profiles (similar in that to real fuel distribution in shaped with RT-instability fiame 
zone in simulations of SN la), which are furthermore insensitive to A. For this model 
we present simple fits of the the parameters entering (1.3) as functions of A, yielding 
very simple implementation of the new model proposed in flame capturing numerical 
scheme. We conclude Chap. 2 by presenting results of actual implementation in FC 
technique, observed fiame speeds, widths and profiles. Discrepancy with steady-state 
results are clearly connected to discretization effects by varying resolution, number 
of cells within the flame width. 

Chap. 3 is devoted to studying noises produced by FC scheme in ID simulations. 
Two new quiet models are introduced, that have finite fiame widths and unique 
eigenvalues for fiame speeds. They are calibrated to yield correct "fiame" speeds and 
widths using the methods of Chap. 2. 

In Chap. 4 we study 2D and 3D behavior of model flames. By modeling circu- 
lar flames in 2D we observe that anisotropic propagation speed is a generic feature 
of reaction- diffusion flames, more prominent at higher expansions. This is a purely 
numerical effect, and is cured by increasing flame width. For one of the new quiet 
flame models, sKPP, this effect is especially pronounced: for flame width of about 
8 cells between reaction progress variable values / = 0.01 and / = 0.99 the flame 
propagates along grid axes by > 5% faster than at 45 degrees to them, for expansions 
corresponding to deflagration in WD at relevant densities of 10 g cm~ and below. 
A free parameter of the second of the quiet models of Chap. 3 is used to eliminate 
propagation anisotropy at all relevant expansions. We observe LD instability, partic- 
ularly severe for sKPP model of Chap. 3 at larger expansions. These 2D effects are 
a strong argument to avoid using sKPP model at densities of < 10^ g cm"'^. Second 
of the quiet models, model B, shows significantly larger critical LD wavelength, and 
slower growth of this instability. At all densities above 3 x 10^ g cm~"^ in SN la 
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problem maximal combined 2D asphericity due to both LD and numerical effects 
stays within 1.2% for simulation times exceeding those used in SN la simulations. At 
lower densities, and in 3D simulations, combined asphericity is larger yet tolerable 
for the "flame" width used, about 5.5-6 cells between / — 0.01 and / = 0.99, which 
is thinner than in original model of Khokhlov (1995)^. We present simulations of 
perturbed planar flames in 2D; observed critical wavelengths for perturbations to 
grow agree with theoretical Acr for LD instability (in Markstein approximation of 
infinitely thin flame with curvature-dependent propagation speed). 

We develop quasi-steady technique for estimating Markstein numbers M (quanti- 
fying flame speed dependence on flame curvature for curved multidimensional flames) 
in Chap. 5, and calculate these for models described in the thesis. Acr for LD insta- 
bility depends on M. Comparison with Acr estimated directly in Chap. 4, as well 
as comparisons of exact M found here with M's estimated numerically provides evi- 
dence that Model B is well-understood, and in its observed behavior physical effects 
dominate over numerical ones. 

We conclude that model B proposed here is best suited for use for Flame Cap- 
turing. 



6. The thinner the "flame" is the better for resolving its small scale features, ultimately for 
better accuracy of the simulation. 



CHAPTER 2 
STEADY STATE ANALYSIS OF ID FLAMES 

In the present chapter we decouple equation for reaction progress variable / in (1.3) 
from the rest of hydro dynamic equations, and study its flame-likc solutions in steady 
ID setting. We write a master equation for steady flame profile f{x) in dimensionless 
form, and find dimensionless velocities d and widths w for a few artificial burning 
rates $(/) used in combustion literature, as well as, in Sec. 2.4, for nonconstant diffu- 
sivity K{f) proposed during this work. Accuracy of the numerical scheme proposed 
for finding d and w is tested against analytic results obtained for one model (original 
one, Khokhlov (1995)) and qualitative results for another one, KPP. Restoring nor- 
malization of the master equation (2.2) leads to a simple prescription for choosing 
scale-factors R and K of reaction rate and diffusivity, such that the flame with coef- 
ficients (2.7) will propagate with prescribed speed D and have prescribed width W, 
for any expansion A. These scalings (2.37) are determined based on dimensionless 
d{A), w{A) found in the first sections of this chapter and fitted in Sec. 2.5 for efficient 
numerical implementation. In the last section we check our results in practice, in 
ALLA code. 

Methods for calibrating flame models presented in this chapter are used for dif- 
ferent models in the rest of the thesis; qualitative results for various possible burning 
rates are used for constructing a new flame model (Model B) in the next chapter, 
which we ultimately suggest for use in FC based on properties of the "flames" it pro- 
duces; eigenvalue problem used in this chapter for finding d and getting flame proflles 
is generalized in Chap. 5 for spherical flames in higher-dimensional problems, to study 
effects of curvature on flames. 

2.1 Physical formulation and numerical method 

Here we study time evolution of the reaction progress variable profile /(r), which de- 
termines heat release in hydrodynamic simulations. In ID it is possible to decouple 
equation (1.3) for / from the rest of hydrodynamic equations for any equations of 
state — provided the burning rate can be considered depending on / only. This is 
the case for 1 stage reactions when Le = 1 (when temperature and reactant concen- 
tration distributions are similar, Zeldovich & Frank-Kamenetskii (1938); relevant for 
terrestrial chemical flames), and for FC progress variable /, by construction. Here 
we demonstrate this decoupling, formulate the steady-state problem and describe 
numerical procedure we use for its solution. 
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2.1.1 Boundary problem for flame speed 

All models studied in this chapter are of form (1.3). One particular type of solutions 
of (1.3) coupled to hydrodynamic equations (as well as of real combustion systems) 
in homogeneous medium is steady ID flame front, propagating with constant speed 
into the fuel {Df is defined here as the speed of the flame in the reference frame 
where the fuel rests). Here we study these particular ID solutions, with physical 
quantities (and /) being functions oi x — D jt only. In such a steady solution it is 
convenient to study the system in a flame rest-frame, in which all physical quantities 
depend on x only (after Galilean transformation). Matter velocity (in this ID steady 
setting) is determined by continuity equation: 

u{x)p{x) — const = —DjpQ. (2.1) 

This further determines f{x) from (1.3), which after this substitution for u{x) reads 

^ p{x) dx dx^ dx^ ^ ^ 

To close the system one needs to express p{x) via heat release distribution, 
dQ/dx = qdf/dx. For this, in isobaric burning, one determines p{x) from enthalpy 
conservation, 

H{P0,P0) + Qf = H{PO,p), 

which in particular provides po/p{x) as a function of / (depending on the particular 
equation of state). This makes (2.2) a second order equation for / only. 

In certain physically interesting situations, e.g. near the center of a WD (main 
contribution to pressure from ultrarelativistic degenerate electron gas), or for ideal 
gas (terrestrial flames) p oc 1/H at constant pressure^, thus (2.2) assumes the fol- 
lowing form, 

which we use for steady-state flame speed and width estimates. 

A = ^ (2.4) 

Pfuel - Pash 

quantifies how the matter expands as a result of burning; with A defined this way 
eigenvalues found from (2.3) are close to true ones even for situations with Hp shghtly 



1. this generally is the case for "7 equation of state", H = (for ideal gases this is the case 

when specific heat is independent of temperature.) 
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nonconstant across the flame. A = qp^^'^^ ideal gas equation of state, 7 denoting 
adiabatic parameter (7 = Cp/cy for ideal gas). 

For solution f(x) to describe a physically meaningful flame profile it needs to 
satisfy physical boundary conditions, 

/(-oo) = 1 (2.5) 
/(+00) - 0. (2.6) 

These may be satisfied only for certain values oi Dj (a parameter in (2.2)); these 

eigenvalues oi D j arc by definition possible fiamc propagation speeds in ID. Cor- 
responding eigcnfunctions, f{x) satisfying (2.3-2.6), arc steady flame proflles; in 
flame capturing these will determine "flame" thickness and fuel distribution within 
the "flame" — approximately, as long as the "flame" segment propagation may be 
treated as ID steady-state, and up to discrepancies due to numerical (discretization) 
effects. 



2.1.2 Numerical procedure 

First, by representing burning rate and diffusivity as products of constant dimen- 
sionful scale-factors and (/-dependent in general) dimensionless form-factors, 

$(/) = i?<l>o(/), K{f) = KKoif) (2.7) 

we rewrite (2.3) in dimensionless form: 

To accomplish this we have introduced dimensionless flame speed, 

d^Df/y/kR, (2.9) 

and dimensionless coordinate 



X = x^jR/K 

along flame propagation. Eigenvalue d of boundary value problem (2.8) with bound- 
ary conditions following from (2.5-2.6) may be found numerically following the pro- 
cedure described next (see Zhiglo (2007) for more detail). 

As (2.8) is invariant under translations in x it can generically be reduced to a 
first order equation by rewriting it in terms of 

df 
dx 
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considered a function of /: 



|(A„(/)p)-d(l4)+5!^ = (2.10) 



This form is used below for quahtative and numerical analysis of the problem. Cor- 
responding boundary conditions are p(0) = p(l) = 0. Eigenfunctions p{f) are non- 
negative. 

In this section we demonstrate the technique for systems with constant diffusivity: 
the original model of Khokhlov (1995) (/g = 0.3 used there) with step-function rate 

^^(R for/o</<l 

[ otherwise ^ ' 

and KPP (Kolmogorov et al. (1937)) 

$ = i?/(l-/) (forO</<l). (2.12) 

Equation (2.10) is integrated by the fourth order Runge-Kutta algorithm starting 
from / = 1, p = 0.^ We use constant grid spacing (A/ = 10~^ for most of the runs), 
except near zeroes of p (the starting point p{f = 1) = 0, and at most one more), 
where it was refined further. 

Integration is actually started at = 1 — 5, with 5 being the smallest refined 
integration step, as p{f = 1) = appears in the denominator of (2.10). Analytically 
found asymptotic expansion is used as initial value for p{fs) (see Chap. 5 for more 
details in more complicated setting in D > 1 dimensions. Using a = D — 1 = in 
Eq. 5.18, and similar equations below that for p for different fiame models, yields 
asymptotics valid for ID fiames.) 

The d eigenvalue is then found for step-function $o(/) by requiring p\ ^—i = 
0. Namely, Newton- Raphson algorithm (see, e.g. Press et al. (1992)) is applied, 

^^^^dd ^"^^ found simultaneously with p{f), ensuring fast convergence. d{A = oo) 
(found beforehand by solving (2.18)) is used as a seed at each new /q value for the 
first A in a row; for subsequent A's the previous one provides seed value for d; 4- 
13 iterations were enough to get d with 10~^ precision. Results are summarized in 
Sec. 5.1. For the case of KPP burning rate the spectrum of steady flame speeds 
is actually continuous, as for the studied (Kolmogorov et al. (1937)) case with no 
matter expansion. We explain qualitatively this feature of the spectrum in the next 
section. 



2. It is imperative to start from / = 1 for parabolic $o(/) (KPP) as a general solution for 
f{x) near x — > +oo (where / ^ 0) is a superposition of two decaying exponentials, and the faster 
decaying one is lost when integrating dp/df, thus making it impossible to satisfy p\f=i = 0. 
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2.2 Spectrum of flame speeds and flame profiles. 
Qualitative and numerical analysis 

Here we show qualitatively why an eigenproblem with step-function burning rate has 
unique solution for Df but the problem based on KPP has continuous spectrum of 
eigenvalues. We describe asymptotic behavior of flame profiles near / = and 1, 
to better understand physical versus numerical features of ID non-steady simulated 
"flame" in Sec. 2.6, and, with this understanding, to propose new flame models in 
Sec. 3.1. Qualitative consideration is confirmed by numerically solving the eigen- 
problem for d, and by numerical integration of (2.10) for several d's. 

2.2.1 Step-function burning rate 

We start with the model described by (2.8) with step-function burning rate (2.11). 
One can observe that an eigenfunction /(x) is monotonically decreasing, in accord 
with physical expectations (see typical flame profiles in Fig. 2.1). More than that, 
f{x) is convex at / < /o and concave elsewhere. Really, the solution of (2.10) at 
/ < /o is p = df{l + //2A), positive and increasing. It is thus enough to show 
that p is monotonically decreasing in (/o; 1) (then p is automatically positive as 
it goes to at / — > 1 as dictated by the boundary conditions). If p were not 
decreasing there would have existed fr G (/o; 1) such that dp/df{fr) = 0. By (2.10) 
p{fr) = {d{l + d'^p/df^{fr) = d/A > 0, thus p would have been increasing in 

some right neighborhood of f-r, and then (2.10) would require that p grew on entire 
(/o; 1), making p{f = 1) = impossible. 

Any self-similar solution satisfying physical boundary conditions must behave as 
follows (for certain xi)' 

/ = 1 Vx < XI, df/dx (Xl) = 0, < / < 1 at X > XI, lim /(x) = 0. (2.13) 

X 

Really, any solution of (2.10) equal to 1 at any point with non-zero df/dx would 
necessarily exceed 1 nearby. While we could allow systems with profiles / taking 
values outside [0; 1], for our particular model with rate $(/) = at / > 1 any 
solution exceeding / = 1 at some x monotonically goes to +oo to the left of such a 
X, as X ^ (so the boundary conditions are not satisfied). Thus / must either 
become identically 1 on some half-line (— oo;xi] or approach f — 1 asymptotically 
from below as x ^ — oo. Behavior of the solution of (2.10) in the region where 
|1 — /I <C 1 coincides with that of the linearized equation. 



f ^ ci- x/d + C2 exp(-dx). 



(2.14) 
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where 



d = + 1/A) 



(2.15) 



The latter does not remain in the vicinity of 1 as x — oo (in fact is unbounded) 
for any values of integration constants ci 2, unless it is differentiably glued to the 
/ = 1 solution to the left of some xi ■ This completes the proof for the behavior near 
/=!• 

Similar analysis in the region where |/| ^1 yields a general solution of the 
linearized equation 

f = CI+C2 exp{-dx), (2.16) 

which tends to at x ~^ ~ 00 if and only if ci =0. From this it follows that any 
solution of (2.10) such that /(oo) = cannot equal zero at any finite point. It has 
an infinite tail, decaying exponentially. 
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Figure 2.1: Flame profiles [i.e. numerically found eigenfunctions). For each /g (the 
ordinate of the curve intersection with x = 0) 3 pairs of profiles are depicted, for 
1/A = 0.05, 4, and 20, larger 1/A corresponding to the curves intersecting x = ^-^is 
at larger angles. 

Summing up, the desired eigenfunction asymptotically behaves like (2.16) with 
ci = as X ~^ +00, and like (2.14) as x ~^ X1+ (xi = ^w/ R/ K), with ci^2 such 
that 

/(Xi) = 1, df/dx{xi) = 0. (2.17) 

For any fixed xi a solution f^^ satisfying (2.17) is unique, and any solution with 
f{—oo) = 1 is a translation of such. In order for this / to vanish at +cxo it must 
belong to another one-parameter subset of solutions. For this to happen there must 
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be one functional dependence among the parameters in (2.10). As it is shown below, 
this actually leads to a unique value for d for any fixed values of the other parameters 
(A and /g here), for which solutions of the boundary problem exist. 

In the q — {=> p{x) — po — const) case (2.10) actually is piecewise hnear, thus 
the above restrictions immediately yield d{A — oo, /q) as the solution of equation 

fod^ = l-e-^\ (2.18) 

(Cf. (2.26); /(xi) is then expressed in elementary functions as A = oo). One can 
write the solution as expansion in /q asymptotically as d^ — j — ^ /'^° — ^ J/^ — 

JO Jo Jq 

(^■^ + + e-3//o + o ((/oe^/-^o)-^); at /o = 0.3 this yields d 2% smaller 

than the value in Khokhlov (1995) (where in g = approximation f{x) was found 
incorrectly); at A's of interest the difference will be more significant. In the limit 
fo^ 1 d vanishes as = 2(1 - /q) + |(1 - /q)^ + J(l - fof + ... (leading term 
agreeing with an estimate in Zeldovich et al. (1985), p. 266). 

Flame profiles found numerically using the proccdTirc described in the previous 
section arc presented in Fig. 2.1. They demonstrate asymptotic behavior as found 
above, and show how the profile width depends on expansion parameter A. 

2.2.2 KPP burning rate 

It is known (see Zeldovich et al. (1985) for detailed discussion) that the spectrum of 
eigenvalues of KPP problem is continuous, comprising all reals above some di. In 
this section we show that the same holds if one includes the term arising from gas 
expansion, find similar spectrum for a wide range of $'s and verify these conclusions 
numerically. 

One can qualitatively analyze the spectrum for burning rate (2.12) along the same 
lines as for (2.11) described in the previous section. Upon linearization in the region 
where 1 - / < 1 (x ^ -oo) (2.8) yields 



/ = 1 - 5+e-^+>^ - 6_e-H A± = |±y^ + 1, (2.19) 

d — d{l K) as before; thus there is a one parameter subset of physical solutions, 
those behaving asymptotically like (2.19) with bj^ = 0. 

By linearizing (2.8) in the region where / "C 1 one gets asymptotic behavior of a 
solution 

/ ^ b+e-^+^ + b-e-^-^, A± - ^ ± ^^-1. (2.20) 
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There are a priori three different situations for drawing further conclusions: 

1) d < 2: any solution getting to a neighborhood of necessarily becomes nega- 
tive. 

2) d > 2: any b-\-,b- > describe physically acceptable behavior, as does a cer- 
tain subset of 6+ < 0, 6_ > 0. Thus one may conjecture that for all d > 2 there is 
an eigenf unction: it belongs to the described above 1-parameter subset of solutions 
asymptotically approaching 1 as x ^ ^-nd on becoming small at larger x it 
behaves asymptotically hke (2.20), exponentially approaching zero as x — > -|-oo. The 
resulting 1-parameter set of physical solutions corresponds to one of them translated 
arbitrarily along x, i.e. there is a unique flame profile for any d (the term "unique" 
as related to profiles is used below in this sense, i.e. up to translations). 

By analyzing (2.8) one expects the / to decrease monotonically, and it is easy to 
see that a solution cannot asymptotically approach any value in (0; 1) as x — > -|-oo 
(say, by linearizing (2.8) near such a value). Thus a solution will eventually get to 
a neighborhood of / = 0, where it will behave as (2.20); unless it becomes negative 
[i.e. have unphysical b±, 6_ < in part; numerical results below show that this does 
not happen, as well as confirm the form of the fiame profile) it will asymptotically 
go to zero, hence being physical. It is the presence of 2-parameter set of solutions 
decaying to zero at x ^ +oo which accounts for the continuous spectrum of d here, 
in contrast with the situation with step-function burning rate. 

To compare side- by-side one may start with some physical (satisfying /(— oo) = 1) 
/ on the left and check if it can satisfy f{+oo) = as well. For step- function 
there is a unique solution (modulo translations) which goes to zero as x — > -|-oo. 
Other solutions asymptotically approach non-zero values; a solution asymptotically 
approaching some ci < /q reads 

fix > 0) = ci + 2A (^(l + -^^yd{l+^i)x _ 1^ ' (2.21) 
(up to translations in x)- As ci increases from to /o derivative ^ 



/=/o 



mcreases 



monotonically from — /grf + 1 j to 0. If for a given d the (unique up to transla- 



tions) / going to 1 as X ^ "OO has ^ j j ^ [~/o'^(^ + 1)? 0] will go asymp- 
totically to the corresponding ci as x — ^ -l-oo, namely it will be exactly (2.21) where 
/ < /o; if its derivative is more negative, it will be described by (2.21) with negative 
ci until it intersects / = at some finite x (and this is not a solution we are in- 
terested in). Easily established monotonousness properties (how the slope ^ 



/o 



/=/o 



varies with d) prove this way the claim that there is a unique d for which the phys- 
ical solution (according to its behavior near f — 1) goes asymptotically to zero as 
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X — > — oo. For the $(/) (2.12), on the other hand, there are no solutions going 
asymptotically to any constant but as x +C)0, but those going to represent 
a two-parameter set of solutions, and depending on the h± in their asymptote the 



4f 
dx 



may take different values, thus making it possible to match this slope of 

/=/o 

the solution physical near / = 1 for a range of c/'s. (In this case /o denotes some 
intermediate value of /, between and 1.) 

3) d = 2: asymptotic solution (2.20) must be rewritten as fix) = (^ + ^0X)6~2^, 
and there is again a 2D domain of physical (6, Bq), yielding / >0, / — i^Oasx^ +oo. 
In this case a meaningful burning profile exists as well. 




100 



Figure 2.2: KPP flame profiles at different A and d. For each d 3 curves are depicted, 
for 1/A = 0.05, 5, and 20, larger 1/A corresponding to the curves intersecting x = 
axis (where / was set to 0.5) at larger angles, and having larger widths. 

Flame profiles f{x) found numerically for $o(/) = /(I ~ /) ^-^'s shown in Fig. 2.2 
for four values of d. These seem to satisfy the boundary conditions for d > 2, whereas 
the profiles for = 1 (each integrated from p(l) = —df/dx\ /=i = 0) intersect / = 
at finite x with non-zero df/dx- Corresponding integral curves p{f) at the same 
values of d and A are presented in Fig. 2.3. Note that for d = 1 p\ 7^ 0; whereas 
for d >2 it was checked that by refining the grid p{0) became correspondingly closer 
to 0, down to p{0) ~ lO""*^*^ at the bulk grid spacing of 10~^ (at finer grid rounding 
errors in double precision reals dominated). This numerically confirms our qualitative 
conclusions about continuous spectrum [2; oo) for d in KPP model. 

One can observe that the flame width grows fast with d. As in the original 
model (Kolmogorov et al. (1937)) one may conclude that only propagation with 
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the smallest velocity is stable (only asymptotes near / = are essential for the 
argument to hold, and these have the same exponential form regardless of A). Of 
course, if the f{x) at some time corresponds to some eigenfunction above, such a 
profile will propagate with corresponding d. But if one considers a process of setting 
up the steady-state propagation, with initial f{x) corresponding to pure fuel on a 
half-line (or / decaying with x faster than the fastest exponent X-{d = 2) = 1 in 
(2.20)) the resulting self-similar front will be the the eigenfunction with the smallest 
velocity. More generally, by considering the evolution of some superposition \1/ of 
the steady-state profiles found above one concludes that amplitudes of all of them 
but one will (asymptotically) decrease in time in favor of the one with the smallest 
d in the spectrum of \I' (they interact due to nonlinearity of the system). As a 
generic perturbation is a superposition of all eigenfunctions, however small it is it 
will eventually reshuffle the profile into that for smallest d. 
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Figure 2.3: Slope p{f) of the profile of KPP fiame, integrated from / = 1, p = 
at A and d as in Fig. 2.2; larger 1/A correspond to smaller p. Note that at = 1 
p{f = 0) is non-zero, in contrast with the d> 2 curves. 



2.2.3 General burning rate 

Observations above may be extended for more general burning rate function as fol- 
lows: 

i) If $o(/) go6s to some positive constant as / ^ there are no solutions going 
to as X —> +00. From a physical viewpoint a system reacting with finite rate in 
the initial state is unstable and self-similar solutions cannot exist. Problems with 
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$o(0+) = u < do have eigenvalues (discrete spectrum). Corresponding eigenfunc- 
tions are identically zero on the right of some X2-, f = ~^{x ~ X2)^/2 + 0((x — X2)^) 
X ^ X2~- sKPP model studied in the next chapter is of this type. 

ii) For $o(/) — 1^1 + o{f) at / — > non-negative solutions going to as x — > +oo 
exist iS d > I^JJi^ fi > for $o(/) to be positive (what is usually assumed in lit- 
erature. See below for < case). General solution getting to a vicinity of / = 
decays exponentially. Analysis near / = 1<^=^/ = 1 — / -Cl then suggests the 
following for — P + fif + o{f) (below d > 1^ is assumed; f{—oo) — 1 required): 

• P > 0, // > — for all d (the above d > 2^/]! assumed) a unique profile exists 
and 3xi :\/x < xi f{x) = 1. 

• z/>0,/i<0 — a unique profile whenever d = d{l + ;^) > (i.e. for 
\]x\ > /x(l -|- 1/A)^ the spectrum is additionally shrunk); again identically 1 on a 
half-line. For d < 2y^—jl there may exist solutions oscillating around / = 1, glued 
to / = 1 on some {—oo;xi\. These violate < / < 1, yet may be used in FC 
in principle. Besides, f{xi+) < 1 and there may exist situations when / does not 
exceed 1 anywhere. 

• z/ < — no physical solutions. For fl < and d < there may be 
solutions glued to / = 1 on some {—oo;xi], but necessarily f{xi+) > 1. 

• P — 0, p, > — a unique solution exponentially approaching 1 as x — > — oo \/d. 
No physical solutions if /x < 0. 

• P = = 0. If $o(/) = qJ^ + o{f^), f > 1, a unique solution exists. The 
solution of (2.10) may be written as p = Jr(l + ^-pr~^ + C'(/2 + /^)), leading 

to / ~ ((1 — f)qx/d)'^~'^ (? > assumed). If $o(/) = in some neighborhood of 
/ = 1 bounded solutions exist but / — ^ 1 at a; — > — oo from above. Again, usable in 
principle even though f{x) is not monotonic. 

• $0 ~ Qf^i but f G (0; 1). / vanishes identically on the left of some xi, in its 

_ 2 

right neighborhood / ~ (^g/2(l + f)(l — f)(x — Xl))^- 

iii) $o(/) = + o(/), A* < (=^ ^o(/) < near f — 0). For any d a unique pro- 
file with /(-|-oo) = 0, exponentially approaching. Behavior near / = 1 is described 
as in ii), depending on $o(/ ~^ but because of the unique profile physical near 
/ = the spectrum of d is discrete (or empty). 

iv) For $o(/) = Qf^ ii^ar / = (g > 0) there are no solutions of (2.10) with 
p(-|-0) = when r < 1 - a situation analogous to i). As a further analogy, when 
g < there may be a solution (unique up to translations); for this 3x2 '■ f — ^ 

2 

at a; > X2, f{x a^2~) ~ {{^2 ~ ~ (2(?" + l))j Using r = 

here yields corresponding behavior in i). When r > 1 on the other hand, there 
are multiple solutions going to (the ODE's pecuharity; general solution decays 
OS f {{r — l)qx/d)^^^^~^^), hence one expects the same behavior as in ii) (apart 
from this power-law tail into fuel), depending on the $o shape near f — 1. When 
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$o(/) = in some neighborhood of / = one expects a discrete spectrum of rf, with 
corresponding eigenfunctions behaving near / = 1 according to the $o(/) ^^^^ 1) 
in ii). 

Summing this up, KPP-hke behavior is quite universal for $o(/) going to hn- 

early or faster at / — > 0. $o(/) = near / = leads to discrete spectrum; positive 
$o(/) decreasing slower than linearly (or not going to zero) as / — * leads to absence 
of steady-state solutions. Eigenfunction f{x) has an infinite tail at / = 1 if $ol f^l 
decreases linearly or faster. 

It should be stressed that these conclusions are based on analysis of asymptotic 
behavior of solutions near / = or 1; as with KPP case considered above the fact that 
there is a 2-parameter set of solutions going to zero as / — 0, say, is not enough to 
make claims as to global behavior of a solution physical near 1 — it may still become 
negative or unbounded near 0. For the cases with claimed continuous spectrum of 
eigenvalues the way the second term in (2.10) damps the solution near / = suggests 
that for any fixed A the spectrum contains all reals above some dj^. Some estimates of 
this lowest eigenvalue (for the case without expansion), and references can be found 
in Xin (2000) (1 am grateful to Lenya Ryzhik for pointing out this review to me). 
Numerical studies of several models (with (r, f) = (2, 1); (5, 1);(0.5, 1);(1, 2);(1, 0.5)) 
agree with general claims of previous paragraphs. A new element appears (in contrast 
to KPP) for more general power law decay of the burning rate near / = 1 and / = 0: 
the lower boundary of eigenvalues d/^ may depend on A. Say, for = f'^ {1 — f) it 
was observed that d = 1 was an eigenvalue for 1/A = {0.05; 4; 20}; d = 0.464 was 
an eigenvalue for 1/A — {4; 20} but not 0.05; and d — 0.215 was an eigenvalue for 
1/A = 20 but not 4 or 0.05. Asymptotic behavior of solutions near / = or 1 is the 
same for any positive d, so no estimates for the dj\^ follow from local consideration 
(near ends of integration interval for the /. Unlike the case for KPP model where 
d\ = 2 is pointed out by local analysis.) 

It does not seem feasible to use KPP-like profiles in FC because of the continuous 
spectrum of the velocities, thus long times of relaxing to steady state, and large 
widths with two exponential tails (thus no way to localize the "flame" compactly but 
still have the steep gradient region well resolved on the grid). 

2.3 Step- function: velocities and widths 

2.3.1 Analytic solution 

Analytic solution of (2.8) with constant diffusivity and step-function $o(/) (original 
implementation, Khokhlov (1995)) as found in Zhiglo (2007) is summarized below. 

Translational invariance is fixed by considering a particular solution with f{x — 
0) = /q. In the region x > $o(/) = (no reaction) and solution is a diffusive ex- 
ponential tail (analog of the preheating zone in standard Arrhenius flames) corrected 
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by matter expansion term in (2.8), only significant where / is not too small: 

-1 



/ = 2A 



2A\ 



dx 



1 1 \ 



V2A /o 



-dx 



(2.22) 



Last approximation is for x ^ ^^'^(^ 2A//o). When this region where / 
exponential decay is perturbed by expansion term is narrower than characteristic 
width of exponential decay (l/d) flame width can be reasonably characterized by the 
slope of the profile /(x) at / = /q; there the slope is maximal, thus this way one 
gets an estimate of the width from below: 



Wa 



\df/dx\ 



/=/o 



n -1 



(2.23) 



For comparison, to quantify how long the tail is compared to the rest of the flame, 
we also employ another width estimate. 



x(/ = /-)-x(/ = i) 



1-/- 

with /_ = 0.1. Analytical expression was found for w^; if /o > /- 



(2.24) 



In 



1 + 2A//- 



1 + 2A//0 

This was computed using solution (2.22) at x > and value 



(2.25) 



(2.26) 



for the rightmost point where / = 1, found directly from (2.8) integrated once over 
X- 

On X £ (xi; 0) solution of (2.8) was found as 



//A = -1 



where 



"='a 



(2.28) 



and 1, K are modified Bessel functions; B is an integration constant (real; another 
constant was fixed by requiring df /dx to be continuous at x = 0). 
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This solution has to satisfy boundary conditions following from (2.5-2.6), 



r /o at (JO = (fk/<o 



(2.29) 



thus leading to a transcendental equation defining sought for eigenvalue for d (enter- 
ing through (Tq and ai arguments): 



lo(i 



/o , 1 



A 3ao 



Notation is Iq = l2/3(cro), Ki = Kj^^3(cti), etc. With definitions above this equation 
has a unique solution for d for every /o, A. 




Figure 2.4: Flame velocities. The nine sequences d{N) correspond (top to bottom) 
to /o = 0.1,. . .0.9 with step 0.1. The worst (/o=0.1) fit (2.36) is shown. 



2.3.2 Results for velocities and widths 

Numerical solutions of (2.30) are presented in Fig. 2.4. Figs. 2.5 and 2.6 show Wq 
and for these two models. Flame profiles as in Fig. 2.1 were obtained by direct 
numerical integration of (2.10) with d from Fig. 2.4. Relative difference between d's 
found numerically and analytically {i.e. by solving (2.30)) is less than 5 x 10^"^%. 
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Figure 2.5: Flame width Wa at the same A and /q. Larger widths correspond to 
larger /q. 



It was that large for /q > 0.7; for /q G [0.3; 0.6] the difference was of order 10~^, 
the accuracy required of the d in numerical procedure (accuracy for solving (2.30) 
was set to 4 X 10~^^ and apparently these, "analytic" errors played no role), and 
monotonically increased to 20 x 10~^ with /q further decreasing to 0.1. Errors 
in widths followed similar trends (and "numerical" widths were consistently larger 
than "analytical" ones, whereas velocities were smaller), though there was additional 
contribution from crude estimation (up to the whole grid) from integral curves x{f) 
obtained from p{f ) via further trapezium integration; the discrepancy in 5 was 
less than 8 x 10-^%. 

For the SN la deflagration problem matter expansion is not large, thus it is worth 
trying to treat 1/A as a small parameter. A first-order correction to the solution 
ho = d{l/A = 0)~2 of (2.18) is (Zhiglo (2007)) 

d-' = " = "0(1 + ^ + 0(A-2)) , In = 5/.„ - (2,31) 

When /o is small as well this may be estimated with the aid of expansion from the 
end of Sec. 2.2.1 for /iq: 



5/0 - e-^/^o ( - + 1 - 6/oJ - e-^/^o + _ _ 6 - 6/0 j + 0{e-'/^yf^), 

(2.32) 
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Figure 2.6: Flame width wi{f{),A) = — xl}=o 1 / '^^^ order of /q (from larger 
to smaller Wf, at 1/A = 20) is 0.9, 0.1, 0.8, 0.7, 0.2, 0.6, 0.5, 0.3, 0.4. 



this is usable up to /q ~ 0.3. Error of using (2.31) does not exceed 1% for 1/A < 
1 (but d = h~^l'^ must be used. Expanding d up to 0(1/A) leads to far worse 
agreement.) At 1/A = 3 the error grows to 6.7% for /q = 0.3 (1.8% for /q = 0.2); 
at A > 0.69 (SN la problem, pf^el > 3 x 10^ g cm'^) it is within 2.1% for /q = 0.3, 
0.16% for /o = 0.2. At these /q both expressions for hi give the same agreement. 

2.4 Alternative flames with finite widths 

In this section we present results for flame model (2.2) with non-constant diffusivity, 
^oif) = /'') £ (0; 1), and step-function rate 2.11. It was studied in Zhiglo (2007) 
for a convenient in FC feature of flame profiles being localized, with no tails (as a 
further improvement over flame model due to Khokhlov (1995) studied in the previous 
section, which produces flames with "preheating" exponential tails, and KPP with 
long tails into both fuel and ash). 

Like in Sec. 2.1 we are looking for a solution of eigenproblem (2.10) with p{0) = 

p{l) = (as before p = —df/dx = —\j K/Rdf /dx). The above KQ{f) leads to f{x) 
C^-smoothly monotonically interpolating between / = 1 Vx < XI = —d{l + 1/2A) 
(fixing /(O) = /o; / = 1 - (x - Xl)V2 + d{x - Xlf /'^ + . . . at x -> Xl+) and 

/ = Vx > X2 (/ = ird{x2 - X))^/"(1 + 0(X2 - x)) at X ^ X2-). 

X2 and the total width Wc = X2~ XI i^ay be expressed in elementary functions of 
/o, A and d for rational r (or as incomplete F— function in general). Values r = 1/2 
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and r = 3/4 seem most adequate. Corresponding widths are 



-XI 

r=l/2 



= ^V2Aarctan Jl^ + d(l + ^), (2.33) 



d V 2A V 2A 



w 3 = 2V4a3/4^-1 

c,3 



Arctan ^^^o/^)'^' + In (/o/2A)^/^ - (2/o/A) + l ' 
l-(/o/2A)V2 (/o/2A+ 1)1/2 



+ d{l + ^). (2.34) 

In the last expression the branch of Arctan to be used is Arctan e [0; tt). 
For any r in the hmit A — > cxd (small expansion) 

(as in the r = case o?(/o, r|A) = h''^/^ = do(l + + • • O""^/^; = etc. 
are now some functions of /q and r); as A — > diverges as 

(2A)^ TT /n^"^ 2A 

a simrr I — r a 



Go sin 7rr 2 Gg 1 ~ 

Wc used d = GqA{1 + GiA + o(A)) as in the step-function model with a standard dif- 
fusion term: d{A) dependence is qualitatively very similar. Fits of d{A) are presented 
in the next section. 

For a range of {/q, r} flame proflles deliver what they were expected to originally. 
Their width Wc may consistently serve as a measure of / gradients, and upon coupling 
(2.2) to the hydrodynamic equations one would really get reasonably uniform heat 
release within the flame width. On the contrary, for models with traditional diffusion 
term ty^ was somewhat arbitrary quantity, most of it (for larger /q and KPP to a 
greater degree) might correspond to "tails" of a proflle, and the heat release in FC 
would remain localized, contrary to the intention to more-or-less uniformly distribute 
it over a few cells near the modeled flame front. Some representative flame proflles 
are shown in Fig. 2.7; they are normalized to unit width (that is, reexpressed in terms 
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Figure 2.7: Normalized to unit width flame profiles with /g = 0.6 (left) and /q = 0.2 
(right). The three curves for each r correspond to 1/A G {0.16; 0.6; 3}, lower near 
/ = profiles for larger 1/A. 



of a new X = (% — Xl)/{X2 " Xl)i resulting supp{df / dX) = [0; 1]). Three specific 
combinations (/q; r) are of particular interest for use in FC: 

• r = 1/2. This has advantage that f{x) behavior near / = is the same as 
near / = 1, thus we are unlikely to introduce new problems (compared to the original 
r = scheme), /q = 0.2 is then convenient as the f{x) shape is least sensitive to A 
in its range of immediate interest, 1/A G [0; 2]. 

• r = 3/4, /o = 0.6. For the A's of interest corresponding flame profiles seem 
most symmetric overall with respect to / i-^ 1 — /; this is perhaps a better realization 
of the above idea: as the whole Wc width is to be modeled on a few integration 
cells this approximate global symmetry seems more adequate to consider than the 
symmetry of minute regions near / = and / = 1. 

• r = 3/4, /o = 0.2. Profile gradients are still uniform enough on [0.1; 0.9], they 
drop to zero in a regular manner within reasonable Ax < 0.25wc. At /q = 0.2 the 
profiles seem least sensitive to A, ensuring consistent performance in all regions of 
the star. 

As another possibility it might seem more physical to write the reaction-diffusion 
equation as 

djpf) 

at 

with density appearing under divergence in the third term, 'W {pKV f) . This case was 
studied as well. Normalized to unit total width flame profiles are more sensitive to A 
than the ones corresponding to (1.3); such a model is also less tractable analytically. 
The model with diffusion term of this form and with K{f) = const was also studied 
as a possible alternative to the original one (2.2). Asymptotic behavior of d and 



+ V(p/u)) = V pKKoif)Vf +pmoif), (2.35) 
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the widths at small and large A are the same as presented in the previous section, 
hence the decision to stick to the original (apparently consistently performing) model. 
Physical transport coefficients change after flame passing, making it questionable if 
putting the p under V really makes the model in this paragraph "more physical", 
in view of the artificial K(}{f) dependence chosen vs effective K{f) for turbulent 
burning (the latter are higher in the ash, qualitatively similar to our KqIJ).) 



2.5 Suggestions for implementation 

2.5.1 Fits of the flame speed dependence on expansion 

Analytically found asymptotes for d{A) (Sec. 2.3.2) at small 1/A do not provide 
enough accuracy to be used for flame tracking in outer layers of WD (although the 
errors are within 1% for /q = 0.2 and 1/A < 2.) Next order in 1/A seems sufficient 
at /o < 0.2 yet computations become too involved in r 7^ case (Sec. 2.4). More 
importantly, flame capturing as presented is a general method, thus it is desirable 
to get results with larger range of vahdity in a ready to use form. In this section 
we present flts neatly interpolating between small and large A regions and then 
summarize the procedure for getting R and K for (1.3) in the SN la simulations. 

- + . (2.36) 

^ ^ l + m2/A (l + m4/A)2 ^ ^ 

with mi. . 4 obtained at each /q to minimize the mean square deviation (with weights 
proportional to actual d{A)) was the simplest flt found. The values shown in the 
graph below guarantee 0.2% accuracy for each /q G {0.01; 0.025(0.025)0.975} and 
A G [10""^; 10^] studied (for any A in fact, as comparison of asymptotes shows). For 
/o > 0.3 agreement is signiflcantly better. I have not succeeded to flnd a simple 
expression for mi. .4(/o) providing good agreement for all expansions; this is due 
to delicate interplay between the two terms near 1/A = 0. Complicated form of 
mi 4(/o) is related to different signiflcance of the two asymptotic regions in the flt: 
for /q = 0.01 the d{A) becomes reasonably linear {d/A ^ const) only at 1/A > 300, 
whereas for /g > 0.5 this happens at 1/A > 2. To overcome this a number of other 
possible flts were tested. In short, 3-parameter flts were signiflcantly less accurate 
(though 3% accuracy was achieved, uniformly on A e (0; 00)), flts with more than 4 
parameters did not yield signiflcant improvements; none of accurate flts considered 
admitted simple expressions for its coefficients in terms of /q. 

This is not a major issue as /q is a parameter one can flx at some convenient 
value throughout the simulation, /q = 0.2 seems most suitable for the flnite width 
flames with r = 1/2 or 3/4, as discussed in Sec. 2.4. For constant diffusivity and 
step-function burning rate of Khokhlov (1995), with new coefficients prescription 
described here, we suggest using /q — 0.2, based on the following consideration. 
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Figure 2.8: Fit (2.36) parameters for r = (the first 4 curves according to the legend) 
and r = 3/4 (next 4). Mi ^3 = mi^2,\fk- 



In Khokhlov (1995) /o was fixed at 0.3, which essentially yielded thinnest fiames 
throughout the A range in the problem (A. Khokhlov, private communication. Note 
that this agrees with our findings, cf. Fig. 2.6.) The width of the "fiame" of the 
original model is W ~ 'Wbifo, A)pAx; /3 = 1.5 was used. Ax is a grid size; it changes 
with expansion A proportionally to wi,{fQ,A), the magnitude of this change may 
be seen in Fig. 2.6. Now, the prescription I advocate normalizes width as well as 
the flame speed"^, thus the logical criterion for /g to suggest would be the fastest 
tail (near f — 0) decay in the profile normalized to unit width wj,, for the case of 
traditional FC realization. This translates to smallest l/w^d. This, however, can be 
made arbitrarily small by taking sufficiently small /q; we suggest to stick to /q = 0.2, 
as this yields l/w^d least sensitive to A in (0.4; 00); fiame profiles normalized to unit 
also exhibit fairly low sensitivity to A. At smaller /o at small A (< 0.5) profile 
shapes on / e [0.1; 1] become rather nonlinear, making the choice of Wf, as a measure 
for "width" more questionable. 



3. We can now also easily quantify the error in flame speed achieved in Khokhlov (1995): the 
values of K and R used there result in actual flame speed larger than the prescribed one by a 
factor of (i(/o, A)v^- It is reasonably close to 1 at small 1/A and /o = 0.3: the flame speed that 
original model produces is 7% smaller in the WD center. However it is ~ 1.45 times smaller when 
p = 3 X 10'^, at 0.5C+0.5O composition. 
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Flame speed when matter expansion may be neglected, a solution of (2.18), value 
of which appears in expansion (2.31) for d at small A, can be approximated as 



do = [(l - /ol-4448exp(0.58058(l - f^^) 

with errors of < 0.64% for any /q. For r = 1/2 and r = 3/4 

1/2 



/o 



0.5 



d0= 2(l-/o)/o^ 



r-l 



provides 1.5% accurate fits. 



/o 


r 


mi 


m2 


m3 


m4 


ed, 10-"% 


0.1 




2.9248 


0.081016 


0.23690 


0.32153 


90 


0.2 




1.9588 


0.14574 


0.26929 


0.42824 


66 


0.3 




1.4620 


0.19389 


0.32605 


0.39759 


11 


0.4 




1.1889 


0.24193 


0.30495 


0.39044 


6.1 


0.5 




1.0312 


0.29358 


0.23117 


0.40545 


11 


0.6 




0.90646 


0.34328 


0.15479 


0.42613 


7.2 


0.9 




0.44817 


0.46634 


0.015034 


0.48381 


5.3 


0.2 


3/4 


1.1683 


0.14139 


0.40205 


0.37378 


56 


0.6 


3/4 


0.81846 


0.34834 


0.14334 


0.42779 


9.0 


0.2 


1/2 


1.3969 


0.14439 


0.35064 


0.39886 


52 



Table 2.1: Parameters of fit (2.36) for d{A) for fiame models with burning rate (2.11) 
and constant K {r = effectively, as in the original FC model: first 7 entries) and 
K — (last 3 lines). The last column shows the maximum relative discrepancy 
between the d{A) and its fit (with mi 4 truncated exactly as in the Table) at A e 
[l/3;20],erf=(|Ad|/d) 

max- 

The fit parameters mi 4 for /q = 0.1 . . .0.9, as well as for the {fo,r} values 
suggested in Sec. 2.4 optimized for A e [1/3; 50] are summarized in Tab. 2.1. 



2.5.2 Prescribing normalizations of dijfusivity and reaction rate for 

propagating the "fiame" 

The strategy for using the results of this chapter in FC is as follows. One picks 
the favorite fiame model; this means a pair {fo,r} for the models considered in this 
chapter (based on step- function burning rate, Eq. 2.11, and diffusivity K — f"). 
One finds corresponding dimensionless steady flame speed d and width w, in each 
computational cell. For models studied in this chapter, this is accomplished by taking 



38 



mi 4 from Table 2.1 corresponding to the {fQ,r} used, (2.36) then defines d based 
on local A (2.4), and (2.25) (or (2.33) for r = 1/2, or (2.34) for r = 3/4) determines 
dimensionless width w {wi or Wq)- Next, one defines needed Dj(A) (based on local 
density, gravity, etc; e.g. (1.1), A meaning the grid spacing). 

K^-f--, R^-f/-, (2.37) 
d w d I w 

(where is a desired flame width, say 4A) will then yield scale-factors for the pa- 
rameters (2.7) appearing in the equation governing / evolution, (1.3). This equation 
with so defined coefficients when coupled with standard equations of hydrodynamics 
will lead to a "fiame" with the desired speed D j and width W within the approxi- 
mations adopted in this chapter (steady ID burning, pressure and Hp approximately 
constant across the fiame, discretization effects neglected). 

2.6 Nonstationary numerical tests 

Here we briefiy present numerical verification of the key results of the previous sec- 
tions, which were based on steady-state consideration. The goal is to sec how well 
the prescribed velocity and width were achieved when prescriptions of Sec. 2.5 were 
used. Flame profile sensitivity to the expansion parameter is also studied here in 
realistic non-steady simulations. 

2.6.1 Flame speed and width 

For the simulations presented below ALLA code developed by A. Khokhlov [Khokhlov 
(1998, 2000)] was used. The simulations were performed in ID without gravity, with 
actual WD equation of state, and heat release corresponding to complete 0.5C-I-0.5O 
burning to the NSE composition at a given density. The ID integration domain 
(called "tube" below) was closed at one end (reflecting boundary conditions); 4 cells 
of hot ash were placed at this end for ignition, in hydrostatic equilibrium with cold 
fuel in the rest of the tube. At the other end outflow conditions were imposed. 

Four models were studied: the three proposed at the end of Sec. 2.4, and a 
model with constant diffusivity and /q = 0.2. For each model the simulation was 
run for 4 densities, p = {3 x 10^; 10^; 3 x 10^; 2 x 10^} gcm"^. These yielded ex- 
pansions 1/A = {1.442; 0.7322; 0.3975; 0.1669} respectively (released heat was q = 
{7.399; 5.736; 4.469; 3.382} x 10^^ ergs g'^). The coefficients in Eq. (2.2) were deter- 
mined from Eq. (2.37) with Df set constant 80kms~-'^ for all runs, W — A cells, d 
was defined using (2.36) with relevant coefficients from Table 2.1, exact expressions 
for Wfj (for Model {fo,r} = {0.2,0}) or total width Wc (for the three finite width 
models emphasized in Sec. 2.4) were used as w in Eq. (2.37). 
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Table 2.2 summarizes the results. D is the measured flame speed, determined as 
D = -J— ^"^ash 14/1 r, o characterize the flame width. These widths were determined 

Pfuel at ' ^i^i'^ 

in the following way (following A. Khokhlov): at each timestep from 2001 to 10000 
a number of cells with / between /^^ and fu was found, and then averaged over 
these 8000 steps. Wi was obtained this way while choosing {fd',fu) — (0.1; 0.9), 
W2 with ifsfu) = (0.01; 0.99), and with (UJu) = (0.001; 0.999). Steady 
state proflle had been established by timestep 500; at timestep 10000 the flame was 
still far enough from the open tube end. For all the results in this Table the tube 
length was 620 km, corresponding to 256 integration cells. The last timestep (10,000) 
corresponded to 1.8-2.6 s depending on density. 



p 


D 


Wi 


W2 




D 


Wi 


W2 






(r; 


/o) = 


(0.75; 0.2) 


(r 


;/o) = 


(0.75; 0.6) 


3 X 10^ 


76.4 


2.69 


3.97 


4.96 


73.2 


2.56 


3.92 


5.20 


1 X 10^ 


77.8 


2.70 


4.02 


4.92 


74.2 


2.59 


3.99 


5.24 


3 X 10^ 


78.5 


2.67 


4.04 


4.88 


75.0 


2.61 


4.01 


5.22 


2 X 10^ 


78.9 


2.64 


4.04 


4.70 


75.6 


2.61 


4.01 


5.12 




(r 


;/o) = 


(0.5; 0.2 


) 


( 


r-, fo) -- 


= (0;0.2: 




3 X 10^ 


76.4 


2.55 


3.77 


4.75 


76.7 


3.18 


5.00 


6.42 


1 X 10^ 


78.2 


2.58 


3.85 


4.76 


78.2 


3.22 


5.23 


6.81 


3 X 10^ 


78.9 


2.57 


3.88 


4.75 


78.9 


3.23 


5.38 


7.07 


2 X 10-' 


79.1 


2.51 


3.91 


1.(31 


79.1 


3.22 


5.19 


7.30 



Table 2.2: Flame velocities, D (in km s ), and widths H^i,2,3 (in cells) at different 
densities p (in g cm""^) for 4 flame models 

To see if statistics was sufficient the same simulation was run with 16 times 
longer tube (divided into 4096 cells) for 160000 timesteps for two of the models. 
Flame speed and the three widths agreed within 0.1% (these were averaged over last 
158000 steps in these long runs) with those in Table 2.2. Therefore larger than 4 
cells for the 3 finite- width models (as well as D distinct from 80 kms"-*^) must 
be attributed to discretization: as / changes from to 1 mostly within 4 cells one 
generally would expect errors in estimating gradients via finite differences to affect 
the ffame profile and propagation speed. To quantify this discretization effect the 
model with r = 0.75, /q = 0.2, p = 3 x 10'' g cm^'^ was run in ID domains 512, 
1024 and 2048 cells long, with D = 80 kms~l and W equal to 8, 16 and 32 cells 
respectively (these D and W. as always, define K and R via (2.37), and are the values 
for fiame speed and width one hopes to get with the simulated flame). The results 
are summarized in Table 2.3. These (together with the flrst quartet in Table 2.2) 
illustrate the trend. In part, for the flame width of order 16 cells or wider the 
becomes smaller than the prescribed total flame width, as it should; the difference 
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between the prescribed flame speed D and the actual value also tends to zero when 
the number of zones in the flame increases. The same holds for different densities: say, 
measured flame speed for the same (r = 0.75, /o = 0.2) model at p = 2 x 10^ g cm"'^ 
is {78.8; 79.4; 79.9} km 8"^ for W = {8; 16; 32}A respectively. 



w 


D 


Wi 


W2 


Ws 


8 


77.3 


5.25 


7.48 


8.40 


16 


79.1 


10.57 


14.67 


15.87 


32 


80.0 


21.32 


29.27 


31.53 



Table 2.3: Flame velocities and widths at different prescribed flame widths W — 
^l/tl- (™ cell sizes. A), (r; /o) = (0.75; 0.2), p = 3 x 10^ g cm-^. 

This study shows that really the difference between the prescribed flame speed 
and the actual one achieved in simulations is due to a small number of cells within 
the flame. The discrepancy can be corrected for by tuning the d and w values in 
(2.37), that is adjusting our analytic prescription for the coefficients in (1.3) or (2.2) 
with additional (A-dependent) factors. This will be done in the next Chapter for 
flame models studied there. 

2. 6. 2 Flame profiles 

Two of the models in Sec. 2.4 were proposed for use for the low sensitivity of their 
flame profiles to the expansion parameter. Fig. 2.9 shows that this nice property 
is not spoiled by the discretization effects, and further clarifies the nature of longer 
profile tails in discretized setting. For each model, /q = 0.2, and r = 0.75 or 0.5, the 
values of / were recorded near the flame position for 4 different timesteps; for each 
timestep the set of values f{xj) was first renormalized in x direction by dividing all the 
Xi by 4A, thus normalizing the numerical profiles to unit total width (more precisely, 
they would have had unit width if there had not been discretization corrections); 
then these renormalized profiles were translated in x direction so that they least 
deviated from the steady-state profiles (Fig. 2.7). This procedure was performed 
with the 2 models displayed in Table 2.2, for the same 4 densities. As Fig. 2.9 
shows the numerical profiles with 4 cells wide flames (1) closely follow corresponding 
steady-state (continuous) profile at all times (apart from the longer tail), and (2) 
as a consequence, are insensitive to density. These numerical profiles, further, show 
about the same density independence for the two models, (r, /q) = (0.5, 0.2) and 
(0.75, 0.2), which suggests to stick to the former one because of its more symmetric 
profile. 

We conclude that all results of this chapter agree for different methods used for 
their derivation: analytical versus numerical solution of steady-state (continuous) 
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Figure 2.9: Theoretical (steady-state; black curve) and numerical flame profiles 
(red dots) normalized to unit theoretical width. Numerical profiles are represented 
by values /(a;^) near the flame position for 4 timesteps at 4 densities. For the 
ij'i fo) = (0.75; 0.2) model the dots are differentiated according to fuel density: circles 
correspond to profiles at p = 2 x 10^ g cm""^, hexagons to 3 x 10^, stars to 1 x 10^, 
and triangles to 3 x 10 . 

problem versus direct numerical simulation using full hydro-solver ALLA. All discrep- 
ancies observed are clearly related to their respective causes: errors of numerical 
integration of ODE's in steady-state approach (effects of which shown to vanish with 
refining resolution; these errors are within 10~^% for flame velocities and widths 
when bulk integration step is A/ = 10~^, further refined up to A/j.gf = 10~^ near 
peculiar points of defining ODE (2.10) — the setup we used for most calculations); 
discretization effects in direct numerical estimates in non-steady setup. Corrections 
due to discretization were seen to be negligible when the "flame" is resolved on 16 
or more cells, but they were up to 10% in obseerved flame speed for some flame 
models when prescribed total flame with W was set to 4 cells, and more than that 
in flame widths. The effect of discreteness of non-steady simulation on flame width 
was seen to be mostly due to longer tails in flame profiles (due to numerical diffu- 
sion). These cross-checks between different approaches provide trustworthy estimates 
for accuracy of corresponding methods, what resolution for numerical integration of 
eigenvalue problem is sufficient, what are discretization corrections for models with 
different flame profiles. The same checks are performed in the following chapters 
for 2 new flame models, for which the methods described in this chapter are used 
for similar calibration, needed for use in FC and for head-to-head comparisons of 
numerical and physical non-steady effects between the models (for which parameters 
like D and W are fixed the same for all comparisons performed). We study numerical 
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noises the flame models introduce in ID simulations in Chap. 3, isotropy of flame 
propagation and flame surface instabilities in 2D and 3D in Chap. 4. 



CHAPTER 3 
FLAMES IN NON-STEADY SIMULATIONS: 
NUMERICAL NOISES IN ID 

This chapter and the next one are devoted to studying features of model flames in 
non-steady simulations. Non-steady properties wc observe (and describe below) are 
visibly separated into physical (thus expected) effects, and effects related to numer- 
ics. The latter are numerical noises, like sound waves generated and flame speed 
fluctuating in time, with periods directly related to flame propagation with respect 
to the grid; anisotropic flame propagation in 2D/3D setting, with grid dictated spe- 
cial directions. The numerical artifacts are clearly undesirable, they disturb and 
systematically deviate numerical model behavior from physically sound one, having 
nothing in common with reality. Finding numerical realization showing least possible 
amount of these numerical artifacts is necessary, and is one of important goals of our 
study. 

Other effects we observe also deviate artificial flame behavior in simulations from 
what would seem ideal for flame capturing applications; effects of this class, however, 
are expectable physical phenomena for any diffusion-reaction flames. These effects 
include dependence of flame propagation speed upon local flame curvature (Mark- 
stein effect, Markstein (1964)); growth of wrinkles on initially smooth flame surface, 
hydrodynamic instability of the type (LD for short) studied in Landau (1944); Dar- 
rieus (1938). These effects depend to a degree on the model diffusivity and reaction 
rate functions; no model, however, is free of these effects. Real physical flame is also 
an example of reaction-diffusion system, albeit more complicated; its propagation 
shows the same physical effects. Ideal flame capturing scheme therefore would be 
not the one showing no effects of this physical class, but the one demonstrating these 
effects with magnitude equal to that of (smeared over grid cell scale) physical flame 
region. Results below. Chap. 4, demonstrate that features of LD instability depend 
strongly on the parameters of the model; we cannot hope that some random flame 
model physical features will match quantitatively those of real physical flame region. 
Detailed study of real flame is required to find its Markstein length, critical length 
and growth rate with respect to hydrodynamic instability; only then one would be 
prepared to analyze deviations from corresponding features of artificial fiame model 
(studied in this work), and try to correct for these deviations. 

Some instabilities, like Rayleigh- Taylor, depend on density contrast across the 
flame zone, and do not depend significantly on specifics of density distribution in 
the transient (flame) region, as long as corresponding instability scale exceeds this 
transient region width. This is not the case for a thick "flame" used in FC, yet 
it is the task of subgrid model to prescribe renormalized effective "turbulent flame 
speed" to correct for this intricate subgrid geometry being smeared out by the thick 
artiflcial "flame". This is not related to the scope of the thesis, this subgrid pre- 
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scription should be the same for any flame model; it may only slightly depend on 
specifics of density distribution within the "flame", which distinguishes one model 
from another (assuming the same widths of the flames they produce). Response of 
the flame to front curvature (Markstein effect) and LD instabihty, on the other hand, 
do depend significantly on specific flame model. They also depend strongly on the 
flame width, thus are expected to be different for the model thickened "flame" and 
for the real thin nuclear (or chemical) flame. What is average manifestation of these 
real hydro dynamic effects on large, resolved scales, whether it is may be compara- 
ble to corresponding effects observed for thickened flame with the width of order of 
characteristic width of the convoluted with instabilities physical flame brush — is an 
interesting topic to study, however this is not touched upon in the thesis. 

With qualitative understanding above in mind, we study several flame models in 
real non-steady simulations without gravity (to have constant fuel density throughout 
the simulation, and to avoid complicating the observations with RT instability, which 
is essentially the same for different flame models). This chapter deals with numerical 
noises which are ID in nature. No physical ID instabilities are known for flame 
models of the type we consider (with only one significant transfer coefficient); we did 
not observe any in simulations performed. 

We describe and cahbrate (using methods of Chap. 2 two new flame models in 
Sec. 3.1. We present results for sound waves produced by these models, as well 
as by the model studied in Chap. 2 (based on step-function burning rate (2.11), 
and power-law diffusivity i^o(/) = = 3/4) in ID simulations in Section 3.2, 

after describing numerical methods used. A feature is identified in the burning rate 
(namely discontinuity), which causes noise in simulations; this noises are too intensive 
for the model used in Khokhlov (1995) as well as in other models studied in Chap. 2, 
based on step-function burning rate. The 2 new models introduced in Sec. 3.1 produce 
acceptable level of noise; this was a criterion for selecting them for further study in 
the next chapters, where 2D and 3D behavior is analyzed. 

3.1 Models studied: definition and steady state parameters 

3.1.1 Definition of the models 

As discussed in Chap. 2 certain features of artificial flame system are desirable for 
flame capturing applications. In part, we want the system to have a unique eigenvalue 
for Dj, and want flame profiles to have finite width (without exponential or power- 
law tails), with reasonably constant df/dx within the "flame". This limits possible 
K{f) and $(/) in (1.3) suitable for flame capturing. 

As a brief overview, two realizations of (1.3) were used in literature before Zhiglo 
(2007). Diffusivity is constant in both, K[f) = K = const, the realizations differ 
by reaction rate forms. The model originally proposed in Khokhlov (1995) has step- 
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function $(/), (2.11). This yields "flame" having an exponential tail into fuel (this 
tail, where / < /q, represents "preheating zone"); D j eigenvalue is unique and exists 
for all R, /g. Another flame model tried in astrophysical literature (Dursi et al. 
(2003)), based on KPP burning rate [Kolmogorov et al. (1937)] $(/) = Rf{l - /), 

has continuous Df spectrum (for any Dj> there exist a solution of (2.3) 

satisfying physical boundary conditions), and two inflnite tails in eigenfunctions f{x): 
f approaches and 1 exponentially as x goes to +oo or —oo respectively. 

In Sec. 2.4 we proposed a new flame model based on the original step-function 
burning rate, but with diffusivity dependent on /: 

m-R^o{f)-[Q tlf<fi^^ ■,K{f)^Kr. (3.1) 

This diffusion-reaction model has a unique steady-state speed D^, like the original 
one (Khokhlov (1995); effectively a limit case of (3.1) with r = 0), and for r > 
corresponding flame profllc is flnite: the region, "flame" , where / is neither nor 1 
has flnite width. We proposed to use this model with 

/O = 0.2, r = 0.75, (3.2) 

as these parameters led to flame profiles having the same shape for any A e [0.25; oo). 
This encompasses all physically interesting range of fuel densities for defiagration in 
a WD, from central density of ~ 2.2 x lO^g cm~"^ (A fa 6.2) down to ~ 3 x lO^g cm""^ 
(considering half-carbon, half-oxygen composition). We will present results for this 
(/o, r) pair only, as this type of model turns out to be quite noisy for any (/q, r) (in- 
cluding original model in Khokhlov (1995), corresponding to pair {fo,r) = (0.3,0)), 
significantly worse than the other two models we will study (and thus of limited in- 
terest for use in flame tracking). We refer to this model as Model A in the following. 

The second model we analyze in this paper is based on a modiflcation of KPP 
burning rate with constant diffusivity, 

$(/) = R{f-ef){l-f + ea) (0 < e/, < 1), 

K{f) = const, ^ ' 

dubbed "shifted KPP" (sKPP for short) at FLASH center meetings where it was 
invented, due to the nature of corrections to KPP burning rate, significant only 
when / is close to either or 1; these corrections effectively cut exponential tails of 
original, "not shifted" KPP-flame proflle (corresponding to ey = = 0) rendering 
flame localized in space, like model (3.1). Besides, this model has unique flame 
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propagation speed, another advantage over original KPP model . We will describe 
the results for this model for ej = ta case: these 2 parameters have similar effect 
on model characteristics due to the symmetry of its burning rate $(/), and there 
was found no significant advantage when using different values in the pair, ^ ea- 
= = 10~^ produces optimal results in terms of 2D/3D behavior and ID noises 
(see below). 

The third model studied. Model B, the one we recommend for use in flame cap- 
turing based on its behavior in multidimensional simulations, is specified by 

$(/) = Rffil-fy^if-ciA)) (3.4) 

K{f) = Kffil-fY-, (3.5) 

with 

= Sa = 0.8, r/ = 1.2, ra = 0.8, (3.6) 

and c(A) e [0.005, 0.3] determined locally as a function of expansion parameter for 
fuel/ash at given local pressure and fuel composition. This model has a unique 
steady-state flame propagation speed (when c(A) > 0); flame profiles are localized, 
and do not have long tails (in contrast with model (3.3)); they look similar to pro- 
files of model (3.1-3.2). Term c(A) was selected to minimize flame propagation 
anisotropy and hydrodynamic flame surface instabilities in 2D at different fuel den- 
sities (Chap. 4). Reasonable results were achieved with c(A) defined as a sphne: 

c = 0.005 for 1/A < 0.515 (3.7) 
c = 0.3 for 0.81 < 1/A < 1.5 (3.8) 
c = 0.2 for 1/A> 1.9, (3.9) 

and continuously linearly changing between these A intervals (sec Appendix for ex- 
plicit formulae). Exponentials (3.6) were chosen based on steady-state considerations 
and non-steady properties. Sf,Sa ^ 0.7 for $(/) to go to fast enough at / — > 
and 1 (to avoid significant ID noises, see Sec. 3.2). Sa < 1 for the model to have 
a unique eigenvalue for flame speed, ry > 1 for the flame not to have infinite tail 
into fuel. Somewhat larger ry, up to 2, produce acceptable results as well, yet no 
improvement can be obtained in 2D behavior. = 0.8 was chosen for reasonably 
symmetric flame profile, with nearly constant df/dx within. Flame profiles for this 
model (Model B), and Model sKPP (3.3) are shown in Fig. 3.4 below. 



1. These properties, as well as asymptotic behavior of the flame profile of this model (namely, 
f{x) is parabolic near both flame boundaries, / = and 1), may be found in qualitative analysis 
in Sec. 2.2.3. 
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3.1.2 Steady-state and numerical calibration of the models: method 

To use flame model in simulations for flame tracking one has to know how to prescribe 
model parameters in such a way so that the model flame propagated with required 
velocity, and had prescribed width. For this calibration we use the same approach 
as in Chap. 2: by factoring burning rate and diffusivity into constant dimensionful 
scale- factors and (/-dependent in general) dimensionless form-factors (Eq. 2.7), 

we get an eigenproblem (Eqs. (2.8), (2.5), (2.6)) for finding dimensionless flame speed 
d = Df/V KR and flame profile /(%) in terms of dimensionless coordinate x = 

x^J r/K along flame propagation. This problem is solved by numerical integration 
of the ODE and using Newton- Raphson method for obtaining the eigenvalue of d, 
for which boundary conditions arc satisfied. The procedure and resolutions used are 
the same as discussed in Chap. 2.1.2. 

To quantify the profile shape, whether the slope of / within the fiame profile 
is approximately constant, or if the profile has long tails, with regions of large and 
small gradients of / having comparable spatial scales, we calculate widths defined in 
3 different ways: 

^'l = x|/=o.9 = X(/ = 0-1) - X(/ = 0-9) 



,0.01 

X 



(3.10) 

Having found d and w (in any convenient sense, e.g. any of wi 2^2) has a simple 
way to normalize the fiame model to yield prescribed fiame speed D and width W 
(same definition of "width" must be used for physical width W for simulations, and 
dimensionless width say. both between / = 0.1 and / = 0.9). Namely (see 
Sec. 2.5.2), one must use (1.3) for fiame tracking, with burning rate and diffusivity 
(2.7), where scale factors are determined through Eq. 2.37: 

Df W Df I W 

K=-f--, R=-f/-- (3.11) 
d w d I w 

To summarize, for any fiame model {i.e. selected for use dimensionless $o(/) 
Kq^J)) one finds corresponding d, and w, and then, using (2.37), finds normalization 
factors for burning rate and diffusivity, R and K, which lead to prescribed "flame" 
(1.3) speed and width. 

Same parameters d and w, defining proper normalization of K{f) and $(/), may 
be found directly numerically, by simulating fiame propagation in ID, finding K and 
R correctly yielding some chosen D and W, and then reversing (3.11) to find d and 
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w; these latter can be used later on for getting K and R yielding any physically 
motivated flame speed. Such numerically found d and w differ slightly from exact 
d and w obtained by solving eigenvalue problem: the former, numerically found d, 
is essentially the eigenvalue of discretized system (2.8), where spatial derivatives of 
/ arc represented by finite differences on a grid, with total flame width of just ~ 4 
grid spacings. Flame profile of discretized problem also deviates from continuous 
one, exact eigenfunction of steady-state boundary value problem. Most notably, one 
gets longer tails near f — and / = 1, thus 2,3 larger than corresponding widths 
of steady-state non-discretized flame (at least for models with well-locahzed flame 
profiles). Effects of discretization on flame proflles and propagation speed for model 
(3.2) were presented in Sec. 2.6; in part, flame speed agrees with steady state one 
up to ~ 1% when total flame width (from / = to 1) was 16 zones or larger; better 
agreement is observed for smaller expansion. 

When flame width is kept constant throughout the simulation (in units of grid 
spacing) it is preferable to use numerically found parameters d and w for discretized 
problem, speciflcally ones found numerically by modeling ID flame with the same 
flame width as the one required for actual multidimensional simulation. This will 
correct for systematic discretization effect (its principal component, which manifested 
itself in ID simulations resulted in d and w deviating from those for continuous, 
steady-state, model). 

For direct numerical (non-steady) estimation of flame speed and width, as well 
as for almost all other numerical studies presented in this and the next chapters, 
dimensionally split piecewise-hnear code ALLA (Khokhlov (1998, 2000)) was used. 
For all runs performed the fuel was taken as half -"^^C, half (by mass) ^^O mixture; 
real degenerate equation of state was used, gravity was set to zero. We assumed 
that the fuel transformed to nuclear statistical equilibrium composition (depending 
on local pressure) and all the heat released in the process was released within the 
"flame". Physically this is not the case at lower densities in a WD (outer layers, 
p ~ 2 X 10 g cm~"^ and below), yet it is the expansion parameter A, dependence of 
flame characteristics on which we are interested in; all the results found are presented 
as functions of A, and not, say, of fuel density (this way, for any flxed A, they must 
be universal for any equation of state, providing assumptions of pressure and Hp 
being constant across the flame are satisfied). For SN la problem at lower densities 
it is thus straightforward to find corrections for our results: one just needs to use real 
heat release within the flame (slightly smaller than what we used, based on burning 
to NSE composition) to flnd correct value of expansion parameter A, and use our 
results (such as d{A), w(A), asphericity of flame surface in 3D, or Markstein number) 
for such found A. 
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3.1.3 Results for flame profiles and velocities 

Results of steady-state calculations of d and w for Models A (3.1-3.2), sKPP (3.3) 
and B (3.4-3.7) are presented in Table 3.1, together with the same quantities found 
numerically. As discussed above, the d and w found analytically, or by finding the 
eigenvalues and corresponding flame profile width through numerical integration of 
(2.2) (we call this "steady-state technique", as opposed to direct numerical simu- 
lations, where physical quantities are discretized on computational grid), will yield 
correct prescribed flame speed D and W (when using (1.3) with parameters deter- 
mined by (2.37) for "flame" evolution) only when W is sufficiently large (~ 16 cells 
or more, as tested for model A in Sec. 2.6). Systematic deviation from prescribed D 
and W grows as W decreases. To correct for this deviation we numerically obtained 
dmim and ^«num yielding correct prescribed Wi and D. These depend on desired Wi 
(tending to steady-state values as Wi — > oo), and should yield the required width and 
speed in 2D and 3D setting (as long as the front may be considered sufficiently flat, 
and providing its corresponding width is close to Wi used to obtain dmxm and Wmira 
in ID), with better accuracy than parameters dgt and Wst found through steady-state 
technique. 

For model A Wi was chosen so that to produce (almost total) width = 6, 
to use the results of Sec. 2.4, where simple analytical expressions were obtained for 
flame speed and total width (/ = to 1) for this model. For sKPP value VFi,nuin — 4 
was chosen based on 2D performance: for smaller Wi surface instabilities are too 
large at densities below ~ 10^ gcm""^. Wi = 3.2 for model B provided reasonable 
performance for all densities studied, from 5 x 10 down to 8 x 10"^ g cm~ . Flame 
speed required in simulations was D = 80 kms"-*^, which is a reasonable value for 
laminar deflagration velocity near the center of a WD (Timmes & Woosley (1992)). 
In outer layers, at smaller densities, laminar deflagration speed decreases, yet one 
typically uses a value of "turbulent flame speed" D r-j St — as prescribed D 

there, which by far exceeds the laminar deflagration speed, due to large gravitational 
acceleration a fair distance from the star center {A denotes Atwood number, Eq. 1.2). 

Dependence of steady state dimensionless speeds d and width wi on expansion 
parameter A is further illustrated in Figs. 3.1-3.2. Shift parameter Sa — ef — 10~^ 
was used for sKPP model. The curves were rescaled to coincide at 1/A = (no 
expansion) to better represent the form of dependence on A. Absolute independent of 
A normalizations are not of much interest for model use in simulations, as prescription 
(2.37) takes care of adjusting K and R for this. The form of the dependence on 1/A, 
on the other hand, is more difficult to account for, and should be fitted on a model- 
by-model basis. Fig. 3.3 shows w^/wi ratio as a function of expansion parameter. 
This ratio is an indicator of flame proflle shapes. This ratio is larger for sKPP model 
(long tails in flame proflle), which is a drawback, as the need to resolve the region 
where gradient of / is large does not allow one to use flame model parameters K and 
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Model 


1/A 


dst 




W3,st 




^num 


^l,num 




^3,nuni 


A 


0.1670 


1.498 


1.225 


1.837 


3.56 


1.424 


1.209 


1.819 


2.035 


A 


0.3980 


1.411 


1.289 


1.919 


3.48 


1.333 


1.268 


1.897 


2.198 


A 


0.7337 


1.306 


1.376 


2.030 


3.44 


1.222 


1.351 


2.000 


2.358 


A 


1.442 


1.140 


1.540 


2.240 


3.41 


1.041 


1.491 


2.184 


2.619 


A 


8.571 


0.5564 


2.401 


3.445 


3.16 


0.4304 


2.014 


2.888 


3.830 


sKPP 


0.1670 


1.681 


8.125 


21.31 


4.0 


1.674 


8.115 


15.99 


21.58 


sKPP 


0.3980 


1.664 


8.893 


23.39 


4.0 


1.596 


8.585 


16.86 


22.82 


sKPP 


0.7337 


1.642 


9.993 


26.44 


4.0 


1.515 


9.319 


18.27 


24.82 


sKPP 


1.442 


1.605 


12.27 


32.83 


4.0 


1.396 


10.91 


21.49 


29.34 


sKPP 


8.571 


1.420 


33.04 


92.19 


4.0 


0.9592 


22.76 


46.77 


65.00 


B 


0.1670 


0.3178 


2.159 


3.424 


3.2 


0.3121 


2.147 


3.441 


4.145 


B 


0.3980 


0.2990 


2.245 


3.595 


3.2 


0.2921 


2.236 


3.583 


4.310 


B 


0.7337 


0.1850 


2.627 


4.096 


3.2 


0.1751 


2.553 


4.058 


4.842 


B 


1.442 


0.1221 


2.865 


4.455 


3.2 


0.1129 


2.744 


4.317 


5.161 


B 


8.571 


0.0661 


3.195 


5.831 


3.2 


0.05816 


3.119 


1.920 


6.016 



Table 3.1: Steady-state flame speeds (dst) and widths (wi,sti '"^3,st); compared to 
the ones found via direct numerical simulations. Wi represents physical flame width 
in number of cells between / = 0.9 and / = 0.1 used in numerical estimates, pf^gj 
corresponding to 1/A shown for a WD with 0.5C+0.5O composition arc 2 x 10^, 
3 X 10^, 10^, 3 X lO"^ and 8 x 10^ g cm~^ respectively. Dnum = 80 km s~-^ for all the 
models. Shift parameters of sKPP model are Ca = ey = 10""^. 
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Figure 3.1: Flame speed scaling as a function of expansion parameter A, found 
through steady-state technique for 3 models, A, B and sKPP with ea = ej = 10"*^. 
The speeds presented were rescaled by their respective values at zero expansion, 
rf^(l/A = 0) = 1.5703, dsKpp{l/A = 0) = 1.6954, rfs(l/A = 0) = 0.33335. Non- 
monotonic dependence of dp on A is due to burning rate dependence on expansion 
through nonmonotonic c(A). 



R yielding too narrow Wi, and the total flame width then, illustrated by value 
is necessarily significantly larger for sKPP model than for models A and B. 

For implementation in simulations using FC, Sec. 2.5.2, we fitted the dependence 
d{A) and wi{A) for Model B. As parameter c(A) in artificial reaction rate for Model 
B is defined as a spline, Eq. (3.7), the fit for speed and width is also given by different 
formulae on different intervals of A; it can be found in Appendix, Eq. (A. 6). 

These fits yield flame speed in ID setting differing by at most 0.8% from required 
80 km/s in a range 1/A G [0.11; 8.6]. Notice that as expansion increases dimensionless 
speed is expected to vanish d{A) oc A, as it was the case for Model A; thus the fit 
should be modified (last part, at large expansions) if it is to be used at 1/A > 8.6. 

3.2 Results for ID numerical noises 

As the flame propagates over the grid one would expect certain amount of acous- 
tic noise generated. Because of discreteness of the grid there is no exact continu- 
ous translational invariance of the problem; say, total burning rate (in simulation), 
summed over grid zones in the flame region (where burning occurs) depends not only 
on continuous flame proflle modeled, but on the precise positioning of that proflle 
with respect to the grid as well. This leads to sound waves generation with period 
corresponding to time needed for the flame to move across one grid spacing. De- 
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Figure 3.2: Steady-state flame width wi = xif = 0.1) — xif = 0-9) scaling as a 
function of expansion parameter A, for 3 models, A, B and sKPP with = ej = 

10""^. These widths are rescaled by their values at zero expansion, wij\^{l/A = 0) = 
1.1762, wisKPpi^/^ = 0) = 7.5687, wisil/A = 0) = 2.0924. 
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Figure 3.3: Ratio of two dimensionless flame widths defined differently, w^/wi {w^ = 
x|j£o^999, wi = xl/=0 9' X is a dimensionless flame coordinate) as a function of 
expansion parameter A, found through steady-state technique for 3 models. A, B 
and sKPP with ea = ef = 10"^. 
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X, cells 

Figure 3.4: Flame profiles for models B and sKPP (black squares), Pf^el = 3 x 
lO^g cm~^, 1/A = 1.442. Numerical profiles were recorded in ID simulations at time 
step 5000. Wigj^pp = 4 and Wip = 3.2 cells, the values used in multidimensional 
simulations. Steady state profiles (continuous line, shown only on the interval where 
/st 7^ or 1.) were obtained by rescaling abscissas of steady profiles by W^i num/'y^l,st 
(so that the rescaled profile width matched numerical one), and then translating in 
X direction for best fit. Note long tails in sKPP profile. 
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pending on boundary conditions these waves may bounce back and forth, perhaps 
seeding and facihtating development of real physical instabilities of the system, thus 
complicating the simulation and making the results questionable. In this section we 
look at local physical quantities distribution in ID domain (simulation tube), as well 
as time dependence of certain integral quantities, like burning rate. 

All models are studied under the same conditions: same computation setup, 
quiet initial conditions (imitating steady state distribution of variables in simulation), 
flame speed is required to equal 80km (which is achieved by choosing flame model 
parameters K and R as described in Sec. 3.1.1. Numerically found values for d and w 
for chosen numerical flame width were used) . Flame widths are different for different 
models, but are required to be equal for each individual model for each run (unless 
specified otherwise) regardless of expansion parameter; these individual widths were 
chosen (for model B and sKPP) based on performance in 2D simulations (see next 
chapter), and are summarized, as W^i^nunn in Table 3.1. 

Flame speed as a function of integration time is presented in Fig. 3.5. It is 
apparent from the plot that while average in time fiame speed is indeed equal to the 
prescribed value for all models, fiame speed for model A fiuctuates by more than 
10% around the average, which is not a desired feature. Also note in this figure 
that transient effects due to not completely perfect initial conditions are completely 
relaxed for all models within 1200 integration timesteps, and that transient deviation 
in flame speed is within 3% for all models. 
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Figure 3.5: Flame speed D{t) = ^p^—^f^ as a function of integration time step. 
Fuel density is 2 x 10^ g cm""^. 
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The pattern of fluctuations is presented in greater detail in Fig. 3.6 for expansion 
parameters typical in SN la problem at the beginning and closer to the end of burning 
in flamelet regime. Relative deviations from prescribed value was rescaled by a factor 
of 200 for model A there for it to have the same scale in the figures as the other 
models. The fluctuations are plotted against physical time, and the range of time 
shown for the two densities is inversely proportional to flame front speed with respect 
to grid, D{l + 1/A). Thus the flame propagates through the same distance within the 
time shown for the two plots, approximately 8 grid spacings. And 8 periods of flame 
speed fluctuations are clearly seen in each figure, pattern of fluctuations in time is 
very similar for each period, but patterns for different models differ drastically. 
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Figure 3.6: Relative flame speed D(t) deviation from its prescribed value Dq = 
80km as a function of time. Fuel density is 2 x lO^g cm~^ (1/A = pfuei/ Pash~ 1 = 
0.1670) for the left figure, and 3 x lO"^ g cm^^ (1/A = 1.442) for the right one. For 
model A this deviation was divided by a factor of 200 so that the rescaled deviation 
had approximately the same scale as deviations for the other two models. 



Next figures show distribution of physical variables in the tube at a fixed moment 
of time (timestep 5000). Matter densities and velocities are presented in Figs. 3.7- 
3.10. Sound waves are visible for all models superimposed on density and velocity 
jump across the flame. Their wavelength is the same for all models at a given density, 
and is different in the fuel and ash (as generation period is the same, but sound speed 
differs). Say, 6 wavelengths are visible in the figures at pf^gi = 3 x lO'^g cm""^ in the 
fuel, between the flame and the right end of the tube. The pattern is complicated 
due to reflections off the boundaries. Reflections off the open (right) boundary were 
observed to be weak, which facilitates relaxation processes. 
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Figure 3.7: Distribution of density and matter velocity in the tube at timestep 5000 
for model A. Fuel density is 3 x lO^g cm""^. Note large fluctuations in matter velocity. 
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Figure 3.8: Density and matter velocity in the tube at timestep 5000 for model sKPP. 
Fuel density is 3 x 10^ g cm""^. 
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Figure 3.9: Density and matter velocity in the tube at timestep 5000 for model B, 
fuel density 3 x 10^ g cm~*^. 
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Figure 3.10: Density and matter velocity in the tube at timestep 5000 for model B, 
fuel density 2 x 10^ g cm""^. 

Distribution of pressure is shown in Fig. 3.11. We show relative deviation from 
average pressure in the fuel: the latter is calculated by averaging over all cells where 
reaction progress variable / < 10~^. Pattern of distribution of the pressure in the 
sound waves resembles that in density and velocity distributions. It is apparent 
that average pressure in the ash (left of the flame) is lower than that in the fuel 
(by ~ /5fuel'"^/c^5 c denoting sound speed). Also note a strong peak in the flame 
zone. Simulations with larger flame width, 20 cells and above, show several short 
pressure waves within the flame, with amplitude rapidly decreasing (as one leaves 
the flame; within a few zones) to the amplitude of pressure waves outside the flame 
zone; for models with step-function burning rate it is obvious from that pattern that 
pressure waves are generated in the cell where burning rate changes discontinuously. 
Heat starts releasing abruptly in the cell when / value reaches /q, generating strong 
pressure pulse. While exact numerical code used, boundary conditions, have effect on 
non-physical, numerical noise behavior in simulations, certain features in flame model 
have universal effect on noisiness of simulations where these models are used. These 
observations on noise generation show that models utilizing discontinuous burning 
rate functions are to produce significant noise. Table 3.2 quantifies our observations. 
On this grounds model A, with 10% fluctuations in flame speed as a function of 
time, and in matter velocity distribution in "steady" ID flame propagation seems 
unacceptable for use in flame capturing. 

Model B has continuous burning rate, which goes to zero as / — >■ or 1. Burning 
rate (3.3) is also continuous, however there is a jump by shift parameter ejoiea when 
/ reaches or 1, which must be a source of noise. Yet with values = Cq = 10^"^ this 
jump is ~ 1000 times smaller than in model A (or original model, Khokhlov (1995)); 
as a result ID noises are comparable for sKPP and model B. We concentrate on 
studying these 2 models further on. 
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Table 3.2: Numerical noise in ID simulations. Each e written needs to be multiplied 
by 10~^ to yield actual relative dispersion. Average flame speed was 80km s"-*^, flame 
width as in Table 3.1 {Wi = 3.2 cells for model B, 4 cells for sKPP, = 6 cells for 
A), except 2 entries for model A, for which W3 = 24 cells. Note that making flame 
wider for model A still does not render noises acceptable. Simulation box was 256 
cells long for expansion parameter 1/A < 1.5 and 2048 cells long for larger expansion. 
For each run dispersion of D{t) was found over N — 8000 steps from 2001 to 10000, 

its relative value Sf) = ^ {D-^{t))t — {D(t))f/ {D(t))t is shown (here (. . .)t stands for 
averaging over the N timesteps). ep represents relative dispersion in pressure in the 
fuel averaged over time: at each timestep relative dispersion of pressure in the fuel 

was calculated, ep{t) — ■\J{p{x)'^)f — {p{x))'f/ {p{x)) /, (■■■)/ means taking average in 
space over fuel cells, defined as cells where progress variable / < 0.0001; these ep{t) 
were than averaged over N = 8000 last timesteps to yield ep shown in the table, ^sh 
represents relative dispersion in matter velocity in the ash in the same manner: at 
every timestep ash(^) = V {u{xf)a - {u{x))l/ {u{x))a was defined ((...) (J meaning 
averaging over cells with ash, where / > 0.9999); these were averaged over the 8000 
timesteps to yield a,sh shown. 
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Figure 3.11: Distribution of relative pressure deviation (from its mean value in the 
fuel po) ill the tube at timestep 5000. Fuel density is 2 x 10^ g cm""^ in the left figure, 
and 3 x 10^ g cm"'^ in the right one. 



CHAPTER 4 
FLAMES IN 2D AND 3D SIMULATIONS 

In addition to numerical noises that manifest themselves in ID as well, there are two 
physical phenomena that should affect fiame propagation in more than one dimen- 
sions. First, normal fiame speed of a curved flame generally depends on the front 
curvature. This is usually called Markstein effect, though this term also refers to 
specific form of that dependence, namely flame speed depending linearly on flame 
curvature: Dji — Dji^^ ^1 — ^f-^ ■ R here denotes radius of curvature of the flame; 

Dji stands for normal propagation speed of a flame with radius of curvature R] pro- 
portionality parameter Ma is called Markstein length. Such linear dependence was 
observed in experiments and simulations (Markstein (1964)). Sec next chapter for 
more detail, and for studying this effect numerically on model flames. 

Second, when fuel and ash densities are different, Landau-Darrieus (LD) insta- 
bility leads to growth of small perturbations on smooth flame surface, resulting in 
wrinkling the flame, and formation of corners later on. Perturbations of any wave- 
length are unstable in idealized case of infinitely thin flame with no Markstein effect 
(Darrieus (1938); Landau (1944)), that is when fiame speed does not depend on 
curvature. When Ma > 0, for infinitely thin fiame, only perturbations with wave- 
length exceeding certain critical value Acr grow. This is what is observed for all our 
models. Besides Markstein length found positive, the fiame width is significant (~ 4 
cells, comparable to initial radius of burnt bubble in simulations). Theory of LD- 
type instability of finite-thickness fiames is much more involved (see, e.g. Matalon 
& Matkowsky (1982)). It can be expected that the critical wavelength exceeds a 
few fiame widths; its exact value, as well as instability growth rate should depend 
on exact distribution of density and heat production within the fiame. Numerical 
effects may additionally distort flame surface. The aim of this section is to study 
how the shape of initially spherical bubble changes with time. 

4.1 Simulation setup 

As for ID simulations of Chap. 3 we use uniform grid with the same grid spacing 
across all our simulations. The simulations were run in square N x N cells box in 
2D, cubic N X N X N box in 3D; ranging from 64 (box side 132 km) to 2048 
(4224 km) for 2D simulations, and from 64 to 256 for 3D ones. Normalizations K 
and R of flame model diffusivity and reaction rate are chosen, as in ID simulations 
described above, via (2.37), to yield "flame" speed 80 km s"-*^ and width Wi as shown 
in Tab. 3.1, unless specified otherwise. 

For all simulations one of two types of boundary conditions is used. In first type 
of simulations defiagration is initiated by a sphere of hot ash placed fully in the in- 
terior of the square/cubic grid (except for a few runs described below initial fiame 
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center coincides with box center); outflow boundary conditions are imposed on all 
boundaries. Initial distribution of reaction variable / and physical quantities is set 
smooth, corresponding to quasi-steady state spherical flame burning outward. Ini- 
tial conditions are characterized by initial spherical "flame" center position {xQ,yQ), 
and radius Rq; the latter signifies here radius of the surface / = 0.9. Explicitly, 
abscissa of ID steady-state flame profile f{x) (steady progress variable as a function 
of dimensionless coordinate, eigenfunction of (2.8)) is rescaled so that flame width 
(/ = 0.1 to / = 0.9) were Wi (the width of the flame we want to achieve) and then 

translated; in precise terms f{x) ^ c(r) : c(r) = / (j,^ ~ ^'o)^^ ^ translation tq 
chosen such that c{Rq) — 0.9. Initial reaction progress variable distribution is then 
set spherically symmetrically around the center of what we loosely call "a sphere of 
burned material at t = 0" , according to this rescaled profile c(r) , r being a distance 
from the flame center. So, if the center of some cell is Rq from the flame center, 
initial value of reaction progress variable is 0.9 in that cell; if it is Rq + Wi off the 
center, progress variable is 0.1 there, etc. / is zero for all cells more than total flame 
width away from the initial "sphere of ash", and is 1 at the center of that sphere 
(unless the flame is too wide, and its tail reaches the center). Other thermodynamical 
variables, and matter velocity are distributed accordingly, mimicking their distribu- 
tion for quasi-steady state spherical flame burning outward from the center (though 
neglecting the difference in flame profile / of spherical and planar flame) . 

Simulations of second type are run in quadrant (2D) or octant (3D): burning 
is initiated from the corner of the cube, by filling a sector of a sphere centered at 
that corner with ash. As for the full cube simulations, this wording means smooth 
initial distribution of thermo- and hydrodynamical variables mimicking quasi-steady 
state 2D/3D burning from the corner, with reaction variable /(r = Rq) — 0.9 and 
/(r — Rq-\- Wi) — 0.1 at i = 0; matter velocity directed radially off the ash sector 
vertex (in computation cube corner), matter resting at the corner. Reflecting bound- 
ary conditions are imposed (throughout the run) on 2 sides/3 faces crossing in that 
corner, and outflow boundary conditions on the rest of simulation box sides/faces. 
Such quadrant/octant simulations are more economical and must provide the same 
results as full cube simulations-*^ , thanks to the problem symmetry about the cube 
center when central ignition is used. These are customarily used in literature for this 
reason. Quite severe numerical artifacts are often observed in such octant simula- 
tions, with system behavior along octant reflecting faces significantly different from 
bulk behavior. We use both types of setup precisely to check how different results one 
gets in full-cube and octant simulations, due to, say, LD instability seeding (or mod- 



1. full cube containing 4 times more cells in 2D, 8 times more in 3D; due to symmetry all 8 full- 
cube octants behave the same, and presumably the same way as one isolated octant with reflecting 
boundary conditions on 3 faces passing through its vertex. 
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ification) by the boundaries, or maybe just imperfectly realized boundary conditions 
— in our simpler setting, with simple chemistry and without gravity. 

4.2 Flame behavior in 2D, theory and observations 

Fig. 4.1 shows how initially spherical flame evolves with time, for the 3 models studied 
with expansion parameter 1/A = 1.442 (pfucl = 3 x lO'^gcm^'^). One can observe 
two features of flame evolution: flrst, some global anisotropy develops, it looks like 
flame propagates faster along axes than along the diagonal of the grid, for Model A 
and sKPP; second, small scale wrinkling of the flame surface is observed, resembhng 
LD instability development. It is seen that Model sKPP, despite having largest width 
Wi (and by far largest 1^3) demonstrates the largest distortion of spherical surface 
as deflagration progresses. The fact that Model A behaves better shows that these 
2D distortions are unrelated to ID noises, by far strongest for Model A. This rules 
out possible explanation of 2D behavior by relating it to flame speed fluctuations in 
time, different in different directions. Pressure waves due to numerical noise are by 
far most intensive for Model A in 2D/3D as well, yet these waves are also not the 
most important factor in disturbing flame surface, as comparison of flame surfaces 
for the models demonstrates. 

These observed phenomena are detrimental for use in flame capturing. Of course, 
real flame region is subject to its own instabilities, it interacts with real physical 
sound waves and turbulence; all of these distort its surface, initially spherically flame 
will quickly become non-spherical in real world, as our model flames do. However, the 
fact that the amount of distortion strongly depends on the model, even for compara- 
ble "flame" widths (making any flame thinner aggravates asphericity development, 
so when sKPP model width is made to match that of Model B its behavior becomes 
further worse) shows that there is no hope that some random model instability devel- 
opment would miraculously match that of real physical flame region. Anisotropy in 
flame speed is clearly a numerical artifact and is not observed in nature. There is no 
control over small scale wrinkles developing on "flame" surface (apart from making 
the "flame" wider; see results below). Therefore the best model is clearly the one 
preserving spherical flame surface as closely as possible, and this is how model B was 
designed. 

Evidence below suggests that small-scale perturbations have the nature of LD- 
instability. It is then understandable that actual distribution of matter density within 
the "flame" in part, influences smallest unstable scale Ad-, as well as instability 
growing rate, which naturally explains why different models behave differently. From 
the figures below it follows that sKPP has the smallest Acr, and instabilities grow 
fastest. Models A and B yield fiames having density distribution within the fiame 
closer to linear, this seems to make them more stable to LD-like instability. As 
observed in Khokhlov (1995) approximately linear density distribution is also the 
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Figure 4.1: Flame for the 3 models, represented by surfaces where reaction variable 
/ = 0.5. Coordinates are in units of box size, e.g. x = ±0.5 corresponds to right/left 
boundary. Fuel density is 3 x lO^gcm""^. The outermost surface is that of Model 
B, at timestep 14000, when flame radius is R — 22Rq; initial radius Rq — 30 km 
for all runs shown in this Figure. This flame is also shown at earlier time, when 
its radius is 9Rq. Innermost circle represents flame surfaces at t = 0. For Model 
sKPP surfaces are shown (thick solid lines) for timcsteps when R = 7Rq and when 
R — 18Rq. Note that global asphericity due to anisotropic flame speed is the most 
prominent feature at early times, and at later times small-scale LD-type instabilities 
severely distort flame surface. Dotted line shows flame position for Model A when 
its radius R — 20 Rq. 
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case for real flame zone (averaged over scales corresponding to typical numerical grid 
sizes) in RT-driven deflagration simulations, in SN la simulations in part. Hence 
there is more hope that Model B (and not sKPP) would behave more similar to real 
flame region in terms of flame perturbations evolution. 

Fig. 4.2 shows flame surface for model B in simulations in quadrant 1024 x 1024 
cells at late times, when radius is large enough to exceed minimal unstable perturba- 
tion length, and enough time has passed for the instability to develop. In this figure 
small-scale surface features evolution, seen in Fig. 4.1 for Model sKPP, looks hke 
typical LD instability development. Notice that, as expected, characteristic length 
of perturbations increases with fuel density, as expansion across the flame decreases. 
In Fig. 4.1, at 1/A = 1.442 just first signs of instability development can be seen, 
with length of order of a quarter of the flame circumference, ~ 400 cells. Thus, even 
though at any density flame with large enough radius will become LD unstable, for 
densities pf^ej ^ 3 x 10^ g cm~"^ time required for instability development exceeds 
typical multidimensional simulation durations. More careful quantitative analysis 
shows (some numbers are presented below) that there is not much difference be- 
tween quadrant and full-square runs (ones with central ignition) in terms of surface 
perturbations as long as flame arc length in quadrant does not exceed a few Acr, 
characteristic lengths of LD-type perturbations. After that boundaries do contribute 
to instability growth in their vicinity, yet for the durations our runs had (enough for 
the flame to pass ~ 1000 cells) there was not seen apparent difference far enough 
from the axes for full-square and quadrant simulations (for Model B). This may be 
seen for a few runs shown in Fig. 4.2. 

At smaller expansions all models behave better in 2D, no asphericity is seen by 
naked eye in simulations with fuel density of above 3 x 10^ g cm^'^ (1/-^ < 0-4; 
some quantitative results are presented in the next subsection). At larger expan- 
sions, however, sKPP flames behavior deteriorates rapidly, making sKPP model use 
in FC dangerous at physically interesting densities (down to a few times 10^ g cm~"^, 
although there is less time for instability development at smaller densities, as corre- 
sponding regions are reached later in simulations with ignition close to star center). 
When sKPP model shift parameters are increased the flame becomes more 

stable to LD-type instability (being another piece of evidence suggesting that long 
sKPP profile tails may be related to predisposition to LD instabihty for sKPP model: 
these tails are shorter for larger parameters e^, e^-); global fiame anisotropy is not 
improved significantly with increasing flame parameters. ID noises become objec- 
tionable at larger shifts, reaching, e.g. 1.8% average dispersion in flame speed at 
ea = ej = 0.1, 1/A = 1.442 (0.54% at = ej = 0.03, 0.38% at = ej = 0.01). 
Values Sa = €f = 10^^ yield ID noises close to those for Model B (at corresponding 
flame widths we use); sKPP is studied with these e^j values below, unless indi- 
cated otherwise. There might be a hope that making sKPP flames slightly wider 
could drastically improve their stability. As Fig. 4.3 shows, LD-like instability is 
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really suppressed for wider flames, expectedly; however, global flame propagation 
anisotropy decreases rather slowly with increasing flame width. Hence we are forced 
to conclude that for use at larger expansions, 1/A ~ 0.6 and larger. Model B is the 
only viable candidate of all the models we considered. 

4.3 Numerical results for flame shape evolution, 2D 

Here we trace how flame shape changes with time for models B and sKPP, for a range 
of fuel densities. Besides comparing model to model in similar conditions we also 
want to compare flame asphericity development for different initial flame radius, and 
for different distances to the boundaries. By varying initial radius we want to better 
understand the role of initial conditions, whether instabilities readily grow from the 
very beginning of the computation (and then this growth will be further accelerated 
in astrophysical simulations with gravity), or whether there is some initial interval 
of time when flame radius is small and LD-type instabilities decay (then, in part, 
unavoidable imperfections of initial conditions will not play that signiflcant role in 
flame behavior at later times, which is desirable); what the dependence of global 
flame propagation anisotropy on flame radius is. To check the role of boundaries on 
flame surface evolution we compare behavior of the models at different resolutions; 
while keeping grid spacing the same this effectively moves boundaries further away 
from flame surface, whereas flame width and initial radius remain flxed, both in 
terms of number of cells and in physical units. 

In Table 4.1 flame geometry is characterized by 2 numbers (called flame aspheric- 
ity parameters below). First, by interpolation we flnd a set of points (r^, (pi) (in polar 
coordinates) on the surface / = 0.5. Among these we flnd the two with extremal 
distances from flame center rniax and rj^in) average radius over ah points, R — (r^), 
and the best fit of distribution with = rQ + ri cos40j. The two numbers we use 
to describe flame surface are Ar/R = (rmax — i^minj/R and ri/rg (AH values written 
in Tab. 4.1 are multiphed by a common factor of 10^.) These numbers give rough 
idea about flame surface: if it is circular both numbers are zero; if ri/rg is close to 
one half of Ar / R one would conclude that large-scale anisotropy of flame propaga- 
tion speed is the main feature, dominating over amplitudes of small-scale LD-type 
wrinkling; if ri / tq is signiflcantly less than one half of Ar / R the wrinkles yield more 
contribution into Ar/R than systematic large-scale anisotropy of flame surface. 

As observed in the previous section, flame shape in flgures, for Model sKPP 
ri/rg > 0, flame surface is in average closer to its center along diagonal of the 
grid than along grid axes. For Model B this systematic anisotropy is controlled by 
parameter c(A) in the burning rate. If, at a given expansion, term c(A) is used that is 
larger than its optimal value (which we approximate as (3.7)) the flame propagation 
speed along grid diagonal exceeds that along grid axes, and vice versa; this anisotropy 
is related to dependence of density distribution within the "flame" on c(A) (governing 
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Figure 4.2: Flame surface for Model B at late times, for high expansion, correspond- 
ing to fuel densities 8 x 10^, 5 x 10^ and 10^ gcm^"^. For each density two flame 
positions are shown (indicated by flame radius written next to corresponding sur- 
face, in km): when perturbations by LD-type instabilities start to be visible, and 
when they are well developed. Notice that characteristic lengths of perturbations in- 
crease as expansion across the flame decreases; also that the flame retains reasonably 
spherical shape (globally) until late enough times: R = 540 km is the same as the 
radius of sKPP flame in Fig. 4.1. Initial radii of the flames were Rq = 120 km for 
1/A = 2.33, Rq = 90 km for 1/A = 3.18, Rq = 30 km for 1/A = 8.57. For 1/A = 8.57 
results of one more simulation are shown (one flame surface when its radius reached 
900 km, short dash). This simulation was run in full square (with central ignition; 
only one quadrant is shown, x = 0.5 corresponds to simulation box boundary for 
this last surface, whereas a; = 1 is right box boundary for rest of the runs), initial 
radius Rq = 240 km. It is apparent from comparison that characteristic lengthscales 
of surface perturbations are determined solely by expansion, and do not show any 
significant correlation with initial radius. No large scale structure anisotropy of flame 
surface is visible for the central ignition run, even though its surface shown is just 
15% of the box size off the box boundary. Finally, at times shown there is not much 
difference between large scale anisotropy for central and corner ignition; it is more 
noticeable at later times, and perturbations generation on reflecting boundaries can 
be seen, but overall the difference between quadrant and full-square runs is not large 
(for Model B) until so late times when open boundaries start playing a role (this was 
checked for full-square simulations in up to 2048 x 2048 boxes for some densities). 
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Figure 4.3: Flame surface for Model sKPP with Wi — A, 8 and 16 cells, pf^gi = 
3 X 10^ gcm""^. One quadrant of full-square simulation (1024^) is shown; x = 0.5 
represents right boundary of the box. Solid line represents / = 0.5 curve in a Wi = 4 
flame when its average radius is 600 km. Dotted line shows the curve defined the 
same way for Wi — 8 flame when its radius reached 630 km. Dashed line shows 
Wi — 16 flame with radius 660 km. Note that LD-type wrinkles are fully developed 
for Wi — 4 flame, but for Wi = 16 flame just first signs of instability development 
are seen at the radius shown; lengthscale of the wrinkles increases with fiame width. 
Also note that the flame is still noticeably closer to the center along grid diagonal 
than along the axes even for Wi — 16. 
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reaction rate distribution within the "flame"). The values we suggest in (3.7) were 
chosen this way, to optimize flame speed anisotropy. 

It was checked in simulations with different resolutions that boundaries start 
playing a role in disturbing flame surface when the flame passes approximately half 
way from the center to box boundary face. Numbers Ar/R and ri/rg at that point 
start deviating from those observed in the runs with larger box (in terms of number 
of cells; the flames have the same width, and are compared at the same radius for 
smaller and larger boxes), by a factor of order 1.4-2, indicating that actual pertur- 
bation started developing somewhat earlier, and have grown to an extent such that 
deviations in flame asphericity parameters become signiflcant. 

Asphericities summarized in Tab. 2.3 suggest that the role of boundary conditions 
is under control, and seemingly the same for Model B and sKPP, and whether the 
flame is stable or not. Numerical noises in the bulk, far from the flame surface are 
comparable for models B and sKPP, all this evidence suggests that it is not some 
intricate interaction with boundaries which makes sKPP flame speed anisotropic, 
but rather local properties of sKPP flame. Small scale ripples generation may be 
more susceptible to background noises in sKPP model, this partly explaining fast 
growing amplitudes of those. Regardless of whether there is any signiflcant contribu- 
tion to wrinkle growth due to waves reflected from the boundaries. Model sKPP is 
unsuitable for simulations with large expansions of matter across the flame. For clar- 
ity, nonetheless, we performed several runs with off-center ignition to verify whether 
instability growth is mostly governed by local physics near the ffame, or whether 
better boundary conditions may, to extent, improve sKPP flames behavior. In those 
full-square runs initial conditions corresponded to quasi-steady spherical 2D burn- 
ing, with flame center not coinciding with the center of the cube. Results for flame 
centered at points (0; 0.2) and (0.1; 0.1) (in units of box size: x = ±0.5 are two of 
the box boundaries in these units) for expansion 1/A = 1.44 are shown in Tab. 4.1 
together with results for central ignition. It is seen that asphericity parameters do 
not differ drastically for these 3 setups. Same holds true for smaller expansions, when 
asphericities are less, as well as for Model B. Fig. 4.4 shows flame surfaces for the 3 
setups: small scale wrinkles look very similar, regardless of relative distances to the 
boundaries. 

In the next section we present some results for growth or decay of a sinusoidal 
perturbation on a planar 2D flame. Those quantitatively conflrm that the small-scale 
instability we observe on circular flames are really of LD type. This suggests that 
the results should not qualitatively depend on a particular code used for simulations, 
as they are a manifestation of a real physical phenomenon, which should be observed 
with any code. Large scale flame speed anisotropy, as we saw, depends to some 
extent on flame width, it is likely to depend on a particular way derivatives are 
estimated by the code (larger stencils used would have effect of additionally smearing 
the flame, thus stabilizing it to some extent), yet the qualitative result that diagonal 
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Model 1/A 


Ar ri 
B 1 ro 1 


Ar ri 
R 2 ro 2 


Ar ri 
R 5 ro 3 


Ar ri 
R 5 ro 5 


sKPPq, iV = 512 0.167 
sKPPq 0.398 
Bq 0.398 
sKPP 0.734 
B, N = 512 0.734 
Bq, N = 512 0.734 


1.36 0.787 
2.16 -1.23 
2.93 -1.74 
8.57 4.00 
3.12 0.331 

2.37 0.246 


3.42 1.38 
3.07 -1.65 
8.61 -4.80 
19.8 9.87 
5.42 1.60 
5.28 1.49 


5.45 2.76 
3.34 -1.12 
11.8 -5.96 
29.4 11.4 

8.06 1.34 
9.85 -2.33 


8.40 4.19 

I. 58 0.083 

II. 4 -5.59 
49.0 18.7 

11.8 -1.88 


sKPP 1.44 
sKPP, Wi = S 1.44 
sKPP, Wi = 16 1.44 
sKPP, e = 0.01 1.44 
sKPP, e = 0.03 1.44 
sKPP, xo,i/o = 0.1 1.44 
sKPP, I/O = 0.2 1.44 


13.3 4.79 
7.09 2.03 
3.58 -0.988 
6.41 1.62 
3.63 -0.811 
14.6 5.21 
11.6 4.23 


54.6 24.5 

37.5 17.1 
26.1 10.2 

30.0 14.1 

16.0 6.40 

50.6 23.9 

57.1 22.5 


99.2 36.4 
41.9 19.4 

36.5 13.5 

39.6 13.2 
13.5 4.37 
96.4 30.8 

116 32.1 


197 53.6 
60.9 23.8 

47.2 17.0 

79.3 22.2 
34.9 3.63 
204 55.0 
314 33.5 


B 1.44 

B, i? = 90 1.44 
B, i? = 180 1.44 
B, R = 15 1.44 
Bq, N = 512 1.44 


8.17 -4.58 
6.86 -3.64 
4.03 -1.78 
7.61 3.69 
7.93 -4.68 


2.77 0.59 

9.47 -4.51 
7.89 -3.61 
17.6 9.43 
2.66 1.09 


3.75 -1.17 

12.1 -5.68 
9.69 -4.41 
8.78 3.79 
2.93 -1.08 


8.80 0.046 
11.0 -3.11 
11.9 4.31 
9.57 3.24 
7.89 0.14 


B, i? = 90, N^2^^ 2.33 
B, = 90 2.33 
B 2.33 
B, i? = 90 3.18 
B 8.57 


12.3 -6.23 

12.4 -6.23 

20.0 -11.7 

14.1 -6.76 
4.31 -2.22 


18.2 -9.20 

17.9 -8.93 
15.8 -7.09 

20.3 -10.6 
6.68 -2.95 


23.5 -9.35 

26.5 -9.25 
18.5 -3.92 
24.7 -11.9 
12.7 -5.05 


50.5 -10.5 

48.6 -1.86 
34.9 0.09 
43.0 -12.3 
34.6 -5.69 



Table 4.1: Asphericity of flame surface at different times, for models B and sKPP. 
Each ^ shown needs to be multiplied with a factor 10""^ to yield actual same 
for ri/rQ. Each pair of asphericity parameters is shown for 4 different moments in 
time, when flame surface radius reached Ri — i?o + 30(km) (the 1^* pair of columns), 
i?2 = Rq+90 (the 2'^*^ pair), R^ = i?o+270 (the 3"^^ pair), and R^ = i?o+600 (the last 
one). Unless specifled otherwise each model is assumed to yield D = 80 km s^"*^, and 
width Wi=3.2 cells for Model B, and 4 cells for sKPP. Initial flame radius Rq = 30km 
(unless different Rq is specifled; Rq and R are in kilometers in this Table.) Default 
setup is central ignition, 1024^ grid cells; simulations in quadrant with corner ignition 
are indicated by letter "q" at model name; if the side of a simulation box is different 
from N = 1024 grid cells, actual is specifled next to model name. For default 
Rq = 30 km the moments of time shownroughly represent moments when boundary 
conditions start playing a role for central ignition simulations in boxes with sides 128, 
256, 512 and 1024 cells respectively (this corresponds to the time when flame surface 
is about its radius away from the boundaries; this is true regardless of whether the 
flame is LD-unstable or not (at radii shown), as may be seen on the results in the 
Table for Model sKPP at 1/A > 0.734 or Model B at 1/A > 2.33. 
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Figure 4.4: Flame surface for Model sKPP with different locations of a seed flame at 
t = 0. pfucl = 3 X f 0'' g cm^^, D = 80kms~^, Wi = 4 cells, Rq = 30 km, each flame 
is shown when its radius reached 360 km (at timestep about 7500). Dashed line shows 
the flame with central ignition, xq — uq — (box boundaries are at x,y — ±0.5); 
solid line corresponds to the flame with xq = j/q = 0.1, and dotted one to the flame 
with XQ = 0,yQ = 0.2. Notice that the first signs of large-scale distortions are seen on 
the fiamc closest to the boundary, due to matter speed asymmetry, yet small scale 
perturbations on all three flames show essentially the same amplitude. 
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propagation is slower than that along the axes (for sKPP) is likely to hold true 
for any code. It would be interesting to systematically check our results with a 
different numerical scheme. While we have not performed systematic code-to-code 
comparison, typical results obtained with code ALLA were confirmed with FLASH code 
(Fryxell et al. (2000)), for Model sKPP developing global asphericity and LD-type 
instabilities with short unstable length, and Model B showing much better behavior in 
this respects at lower fuel density. Derivatives are estimated with 4*^ order accuracy 
in FLASH code (vs 2"^"^ order in ALLA). Fig. 4.5 shows a typical flame surface simulated 
in quadrant 256 x 256 with FLASH code, with the same flame model parameters as 
discussed above, with 1/A = 1.44. Same features of the model are observed. 

,|: I I I I I I '3 




1vlo' 2vlo' JiflO'' 4vlo' 51^10' 6->^10'' 

Figure 4.5: Flame surface for Model sKPP, Pf^el = 3 x 10^ g cm""^, D = 80 kms 
Wi = 4 cells, computed in quadrant 256 x 256 cells with FLASH code. 



4.4 Perturbed planar 2D flames, LD instability 

Here we present results of a study of the instability development on sinusoidally- 
perturbed planar flames in 2D. We simulate such flames in rectangular domains with 
reflecting boundary conditions on 3 of the sides, and outflow boundary condition on 
the fourth side, towards which the flame propagates. We observe that for domains 
with a width below certain critical value Lcr (which depends on expansion parameter 
and the flame model used) the perturbation decays as the flame propagates. For wider 
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domains certain non-planar flame surface develops, which steadily propagates into 
the fuel. 



Theoretical background 

The behavior we observe, flames in wide domains developing non-planar front shape, 
qualitatively agrees with known from experiments behavior of terrestrial flames; how- 
ever, no analytical studies have been done in the regime relevant in FC, with fairly 
wide flames and non-negligible expansion of matter across the flame. 

First theoretical studies of this type of instability, linear analysis of its develop- 
ment (when the perturbation is small), were done by Landau (1944) and Darrieus 
(1938) in approximation of infinitely thin fiame, with flame velocity independent of 
its geometry and background hydrodynamical flows; the result is that perturbations 
with any wave number k grow exponentially, with the rate cu proportional to k: 

^Ld(^) = kD — . (4.1) 

cr -I- 1 

D here denotes flame speed, a = Pfuel/Pash = 1/A + 1. As discussed above this 
approximation is invalid at small perturbation lengths, comparable with the flame 
width. Real flames are generally stable with respect to such short-wavelength per- 
turbations. 

A simple (and physically sound) modification of the original LD consideration 
leading to stable perturbations with large enough k is to take into account fiame 
speed dependence on the fiame curvature. Assuming linear dependence on curvature 

1/R, D{R) = Dji=(yQ ^1 — "^j, stfll in the limit of infinitely thin fiame, modified 

dispersion relation is given by (Markstein (1964)) 



kD 



(7+1 



-cr(l + Ma A;) + \J a"^ + a"^ - a + {Ma k - 2a)Ma ka'^ 



(4.2) 



For positive Ma, which is the case for all fiame models we consider for FC applications 
and for most chemical fiames in practice^, perturbations with 

27r/fc < Acr(Ma) = 47r(A + 1) Ma (4.3) 

decay. This is understood qualitatively, as positive Ma stabilizes the fiame by de- 
creasing propagation velocity on convex parts of the fiame surface (the ones that are 



2. Ma may be negative in, e.g., systems with significantly different diflusivities of different 
reactants when the most diffusive reactant is deficient (for example for lean hydrogen-oxygen flames), 
cf. Matalon et al. (2003) for theoretical study. Infinitely thin fiames with negative Ma are unstable 
for any perturbation wavelength, as are flames in original LD approximation, with Ma = 0. 
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displaced forward into the fuel, with respect to a plane unperturbed flame surface), 
and increasing the propagation velocity into fuel on concave parts, which are behind 
the unperturbed surface. This effect tends to suppress growth of the instabilities 
(due to main component of LD instability, with D{R) dependence neglected), and 
it does stabilize the smooth flame surface at sufficiently short perturbations (when 
curvature is larger, and this curvature effect is more pronounced). 

Linear analysis was also performed up to the flrst order in flame width (in Mat- 
alon & Matkowsky (1982); Pelce & Clavin (1982); Matalon et al. (2003), under 
different conditions), thus elucidating the effect of flames being not infinitely thin, 
and obtaining some estimates for Markstein lengths. These works only deal with 
Arrhenius-type reaction rates, thus the results are not directly applicable to our flame 
models; besides, especially for sKPP, the "flame" width is quite close to observed 
minimal unstable wavelength, thus the first order terms in asymptotic expansion in 
small Wk are not apriori likely to provide a trustworthy description of instability 
development'^. According to our simulations (summarized below) the critical width 
of the tube, at which stable regime of flame propagation changes from planar front 
to certain nonlinearly stabilized shape, closely corresponds to an estimate based on 
Eq. (4.3) for model B; this quantitative agreement proves that the small-scale insta- 
bility we observed in 2D in the previous section is of the hydrodynamic, LD-type 
nature. 

4.4-2 Numerical results, comparison to theoretical estimates 

Results for critical tube side Lcr orthogonal to flame propagation, its smallest value at 
which initial perturbation develops into certain nonplanar shape (instead of relaxing 
to planar surface, which happens at smaller L), arc summarized in Tab. 4.2. We 
used 512 cells long computation domains for larger expansions, corresponding to fuel 
densities in a 0.5C-I-0.5O WD of 10^ g cm~"^ and below, and 2048 cells long domains 
for smaller expansions; initial perturbation was one sine wave on the domain width 
for most of the runs, with amplitude equal to O.IL. Initial distribution (in the flame 
propagation direction) of progress variable / and physical quantities corresponded 
to ID steady flame proflles. For larger expansion, 1/A > 1.44, uncertainty in Lcr 
is about 2 cells, the smallest amount by which we could change the tube side; the 
smallest supercritical width observed is shown in the table. At these larger expansions 



3. Markstein numbers, which are the ratios of Markstein lengths and corresponding flame widths, 
are typically between —1 and 6 for chemical flames, according to the cited theoretical estimates, 
and to experimentally found valuers. Markstein numbers for our model flames arc between 0.2 and 
0.5 for Model B (larger values corresponding to larger expansion), and about 3 times smaller for 
Model sKPP with = e/ = 10~^ (and go to zero as sKPP shift parameters Ca, e/ go to zero.) See 
next chapter for more detail. 
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the steady flame shape significantly deviates from planar (see Fig. 4.6 for example), 
and the critical tube width is easy to recognize. 

T/k 0.167 0.249 0.398 0.734 1.44 2.33 3.18 3.95 5.05 8.57 

Acr(Ma) 60.3 43.4 30.6 30.8 30.0 19.7 19.0 18.8 18.7 19.3 

L*r 68 48 32 34 32 24 26 20 16 16 

L*,/5, Wi = 16 17.6 16.0 14.4 13.6 

Table 4.2: Minimal tube width (in cells) at which planar fiame is hydrodynamically 
unstable, at different expansions 1/A. Acr(Ma), the 2^'''^ row, is a theoretical predic- 
tion (4.3) in Markstein approximation (of infinitely thin fiame, with linear Markstein 
law for fiame speed dependence on the front curvature); the values of Markstein 
length were taken from Tab. 5.1. The third row shows estimates of the critical tube 
width from above from direct numerical simulations (evolution of initially sinusoidal 
front in 2D tube); these are the smallest values we found to yield non-planar (non- 
hnearly stabihzed) steadily propagating fiames. The estimate from below is 2 cells 
smaller for 1/A > 1.44, and 4 cells smaller for the lower expansions shown. Flame 
speed is 80 kms~^, flame width Wi = 3.2 cells. The last row shows same upper 
estimates for L"^^ obtained for Wi = 16 cells wide flames (V = 80 kms~^), to un- 
derstand whether somewhat different relation between L^j. and estimate Xcv{Ma) at 
higher expansion represents a physical effect (corresponding to increased role of fi- 
nite fiame width), or numerical discretization effect (observed in Sec. 2.6 to be more 
pronounced at larger expansion, for various quantities). This latter L^j. was divided 
by 5, so that to offset the effect of Markstein length Ma (and thus corresponding 
Xcv{Ma)) being 5 times larger for this extra wide "flame". 

For smaller expansion, on the other hand, nonlinearly-stabilized non-planar shape 
at near-critical tube width only deviates by a few cells from the plane, for tube widths 
of tens of cells. This deviation scales approximately as 1/A when this expansion 
parameter is <C 1, in accord with analytical studies (Sivashinsky (1977); Michelson 
& Sivashinsky (1977)), this making recognition of critical tube width less certain. 
Besides, relaxation times become significantly longer at smaller expansion: linear 
growth rate in LD approximation scales as 0.5DA;/A; Eq. (4.2) yields growth rate 
scaling oc l/A"^ (assuming Ma going to some constant value at zero expansion) near 
critical tube width Lcr = Acr(Ma) (4.3): cj^//g(Lcr + 5L) = ^qJ^^2 • The newest 
version of the code allowing stable computation for over 3000C)0 timesteps required 
for densities of 2 x 10^ g cm~"^ and above became unavailable to us in the middle 
of this study, thus uncertainties for Acr at smaller expansions are estimated to be of 
order 4 cells. Only one run was finished at fuel density 5 x 10^ g cm^*^ (V-^ — 0.118, 
XcriMa) = 82 cells); the domain width of 72 cells was found subcritical. 

It is interesting to observe agreement within about 10% between numerically 
found values for the critical tube width and theoretical estimate in Markstein approxi- 
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Figure 4.6: Perturbed planar flame evolution at 1/k — 1.44 for subcritical L — 30A 
and supercritical L = 32A domain width, x and y are in units of the domain length, 
which is 512 cells; y-scale is the same for both L. Flames arc shown at timesteps 0, 
7000 and 14000. The flame at t = 7000 was translated by —0.25 in x-direction, the 
one at t = 14000 was translated by —0.5, to show all three on the same plot. 

mation at expansions (1.44 and below), corresponding to fuel density of 3x lO^g cm"'^ 
in SN la problem. It is even better at expansions 1/A < 0.734 if Markstein numbers 
found in direct numerical simulations (see Tab. 5.1) are used. At higher expansion 
agreement is somewhat worse. The fact that at smaller expansions observed Lcr is 
somewhat larger that theoretical one may be attributed to numerical diffusion tails in 
/ proflles, thus yielding somewhat larger Markstein length than the one found in con- 
tinuous quasi-steady state technique. At smaller expansions, however, numerically 
estimated Lcr becomes smaller than Acr(Afa). To check whether we understand the 
effect of "flame" discreteness properly we reran some simulations at small densities 
with distinct normalization of the "flame" governing parameters, K and i?, to yield 
"flames" with 5 times larger width, for which the effects of discreteness must be less 
pronounced. As Markstein length (and thus Acr(Ma)) scales proportionally to the 
flame width we divided observed upper estimate L^j. for such wider flames by 5, for 
direct comparison with the rest of the numbers in the table. The effect of changing 
the ffame width is signiflcant (as it was for other quantities related to the "flame" 
in Sec. 2.6 at large expansions), yet in the predicted direction: critical perturbation 
length decreases when discretization effects become smaller. Thus we are tempted to 
conclude that the discrepancy observed at higher expansions is due to real physical 
deviation from Markstein approximation, due to flnite flame width (moreover, rep- 
resenting a larger fraction of Lcr at these larger expansions); this deviation seems to 
monotonically increase with expansion. 
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For model sKPP with = = 10^ and Wi = 4 cells the critical perturbation 
wavelength is about 2.5 times smaller at small expansions (in Markstein approxima- 
tion — as are corresponding Markstein lengths, cf. Tab. 5.1 for Markstein numbers). 
These are even smaller at larger expansions, 7 times smaller than those for Model 
B at 1/A = 1.44. We really observe steadily propagating non-planar sKPP flames 
(corresponding to a full wavelength, not reorganizing into a half-wave configuration) 
in domains 10 cells wide at pf^gj = 3 x 10^ g cm""^, 8 cells wide at pf^gj < 5 x 10^; 
Lcr ~ 8 at pf^gi — 10^ — a seemingly stable non-planar configuration in 8 cells wide 
domain starts transforming into a half-wave configuration after running through 120 
cells (this flame relaxes to a stable half-wave configuration when it runs through 250 
more cells). Flame profile, / distribution across the fiame, in such stationary fiames 
is clearly distorted; the fiame is wider where it is convex towards the fuel. At smaller 
expansions Lqt is noticeably larger than according to Markstein approximation, say 
it is 30 cells at pf^gi = 3 x 10^ g cm""^, almost the same value as for Model B (with 
standardly used for the latter Wi — 3.2 though, i.e. thinner than the sKPP flame). 

4.4-3 Discussion 

Steady flame shapes we observe in supercritically wide domains visually resemble 
the ones reported in Rastigejev & Matalon (2006), the only study we are aware 
of dealing theoretically/numerically with non-linear stabilization of flames with fl- 
nite expansion. Results for flame shapes are compared there to analytical results 
of Michelson & Sivashinsky (1977) for small expansion, with convincing agreement. 
Methodology of Rastigejev & Matalon (2006) is quite different from ours: the authors 
try to study non-linear stabilization in Markstein approximation, with inflnitcly thin 
flame (tracked with Level Set technique), Markstein effect introduced by hand in pre- 
scribing flame speed as a function of locally computed curvature. Results for critical 
wavelength are not reported there with reasonable accuracy, although for a few ex- 
pansions it is noted whether existence of stable non-planar shape is in agreement with 
linear stability analysis, Eq. (4.3). Our observations agree with those in Rastigejev 
& Matalon (2006) in that the most stable flame shape has exactly one maximum for 
supercritical tube width; this also agrees with analytical results in small expansion 
case^. 



4. These observations do not seem conclusive though. Say, our results might suggest that a 
configuration with several wavelength corresponding to a stably propagating flame shape, arranged 
end-to-end on a tube width, may be stable with respect to sufficiently small perturbations. Such 
configurations were obtained by starting with several sinusoidal wavelengths on a tube width as 
an initial condition; this evolved to described configuration, rather that to one-maximum one; this 
configuration in several runs propagated (quasi-)stably through entire 2048-cell long domain. More- 
over, certain small perturbations not respecting the discrete symmetry of multi-wave configuration 
seemed to decay, without transforming the multiwave configuration into an absolutely stable one 
with one maximum on tube width. In other runs, however, with superficially similar initial condi- 
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We used reflecting boundary conditions on long boundaries of the computation 
domain, whereas Rastigejev & Matalon (2006) use periodic ones. This makes it pos- 
sible for us to observe another stable propagation mode, different from that shown in 
Fig. 4.6, namely, the one with a half of the steady flame shape occupying the whole 
tube width. Such a mode is expected to be more stable than the one with the whole 
profllc (symmetric with respect to the longitudinal symmetry axis of the tube). How- 
ever the times needed for transition from a (quasi-) stable whole-wavelength shape 
to stable half- wavelength shape, with just numerical noises in our simulations as a 
trigger, were usually long; only in a few runs could we observe such a transition, in 
tubes with width not sufficient to support a stable non-planar profile with one wave- 
length, but wide enough to accommodate a half of the critical wavelength. Since this 
happened more readily at larger expansions, for some of the corresponding results 
we report in Tab. 4.2 we used asymmetric initial conditions, with half of a sine wave 
on tube width (maximum on one side, minimum on the other). 

To sum up, results of this section show that the small-scale instability we observed 
in 2D simulations is indeed a hydrodynamic instability of LD-typc. It is masked for 
Model sKPP with numerical artifacts, like global flame propagation anisotropy, and 
demonstrates significantly different characteristic instability lengthscales and growth 
rates in different directions. For Model B numerical artifacts are almost absent, and 
characteristic for LD instability shapes are easily recognized on 2D flames. 

Dominating lengthscales of perturbations of cylindrical flames show somewhat 
stronger dependence on expansion parameter than observed above critical wave- 
lengths of perturbations on planar flames (we refer to Model B through the end 
of this section). The weak dependence of the latter is in accord with linear results 
in Markstein approximation: Eq. 4.3 suggests rather weak dependence of Acr at ex- 
pansions 1/A > 2.33: dependence through factor 1 + A (weak at small A) is further 
offset with (accidentally?) correlated dependence Ma(A). Theoretical (linear) re- 
sults for LD-type instability for spherical flames (see, e.g., Zeldovich et al. (1985), 
Ch. 6) in Markstein approximation help to qualitatively understand characteristic 
lengthscales of the perturbations, as well as timescales associated with their devel- 



tions (each wave having supercritical length, in part) perturbations started to grow asymmetrically, 
with the shape transforming to a one-wave or a half-wave configuration. Rastigejev & Matalon 
(2006) started with random flame shape, with no symmetries, in simulations dedicated to getting 
one-wave configuration in the process of instability developing towards its nonlinear stabilization; 
thus no results are presented to judge about local stability of multi-wave configurations. This is 
in contrast to analytical results in small-expansion analysis, where instability of any configuration 
distinct from a one-maximum one is proven mathematically. Those results, however, are obtained 
for basically a different problem, with all terms beyond first order in 1/ A dropped (similarly to our 
estimates for d{A) based on expansion in 1/A at the end of Sec. 2.3.2). It is plausible that for the 
full problem there may be islands of stability in perturbation space near multiwave-configurations, 
vanishing as 1/A ^ 0, or perhaps islands existing only when finiteness of the flame width is taken 
into account. This requires further study. 
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opment. In contrast with nonlinear stabilization of flames in the tube, where the 
stable shape (at least in the limit of small expansion) steadily propagating has 1 
maximum, for spherical flames only higher harmonics are unstable; the ones repre- 
sented by spherical functions P^(cos ^)e^"^^ with index n less than certain critical 
value grow slower than the flame expands, and do not lead to the surface distortion 
as time progresses (all harmonics grow according to a power law in 3D and 2D, not 
exponentially as perturbations on a planar flame). This critical index is ncr = 12 for 
positive Markstein length and 1/A 3. For larger and smaller expansion the critical 
index is larger. 

For any supercritical harmonic n a perturbation having corresponding distribu- 
tion over the flame surface will grow after the flame has reached certain radius r{n). 
Markstein effect will not inhibit growth of any supercritical harmonic at late enough 
times, as the linear size of corresponding perturbation scales proportionally to the 
flame radius, and will eventually become large enough for Markstein stabilizing ef- 
fect to become non- substantial. The harmonic that becomes unstable at the smallest 
flame radius will dominate for some time, however we are not aware of any the- 
oretical studies for characteristic perturbations structure at late times. According 
to flame surfaces in Fig. 4.2 for the times studied linear perturbation lengthscales 
stay approximately constant with time, cells expanded with the flame keep subdi- 
viding into smaller cells. Characteristic linear perturbation scale is about 50 cells for 
1/A = 8.57, and larger for smaller expansions, according to the flgure. Ratio of this 
scale to Acr for planar flames is about 3 at this expansion, and somewhat larger at 
smaller expansions. 

To reiterate, theoretical results we summarized correspond to a spherical flame 
in 3D in Markstein approximation, thus spcciflc numbers, like ncr, first harmonic to 
start growing, etc. are not expected to describe well results of our simulations of 
cylindrical flames in the previous section. Qualitative picture, however, presented 
in the previous paragraph seems to apply to 2D as well as 3D setup (following the 
derivations in 3D), and agrees with the features of small-scale instability development 
we observed in numerical simulations. One particularly interesting number to check 
is a critical radius of the flame when the flrst unstable mode appears. According 
to linear analysis in 3D in Zeldovich et al. (1985) this radius is rev = Mar*, where 
dimensionless r* depends on expansion, and is about 60-70 for 1/A e [2.33; 8.57], 
about 80 for 1/A = 1.44, 120 for 1/A = 0.734, 300 for 1/A = 0.398, > 600 for 
1/A = 0.167. For 1/A = 1.44 this yields Tcr ~ 110 cells, or 230 km; this seems 
to reasonably well agree with numerical results in 3D we present in the following 
section. 
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4.5 3D results for flame shape evolution 



LD-type wrinkling of flame surface is expected to grow faster in 3D and may become 
an issue for Model B at larger expansions; it is likely to remain a problem for Model 



Our approach is the same as described for 2D simulations. Flame propagation is 
initiated either from the center of the cubic computation box with outflow boundary 
conditions on all faces, or from the corner, with reflecting boundary conditions on 
the 3 faces passing through that corner, and outflow boundary conditions on the 
other three faces; we refer to the latter setup as a simulation in octant. Initial 
distribution of physical variables and reaction variable / is spherically symmetric 
around the flame center; a sphere (a sector for octant simulations) around the center 
contains hot ash of the reaction at rest; away from the initially burned region physical 
variables correspond to fuel, moving radially outwards, according to quasisteady 
solution for a spherical flame in 3D; in the intermediate region, the "flame" itself, 
physical quantities correspond to a mixture of fuel and ash, / has some intermediate 
value monotonically changing from 1 in the burned sphere to in the fuel ahead of 
the flame, according to ID steady-state proflle. The distance from the initial "flame" 
center to the point where / = 0.9 is called initial flame radius Rq. 

The largest computational box we use for 3D simulations is 256^; LD-type m- 
stabilities do not have enough time to develop to the extent close to that in 2D 
simulations. To characterize flame geometry we use several numbers akin to those 
we used in the previous section. First, Ar/R (deflned the same way as for 2D, by 
dividing the difference Ar between extremal radii on / = 0.5 surface by the average 
distance R from the flame center over all points on the surface) characterizes overall 
sphericity. To get an idea about large scale deviations from the spherical shape, 
and in which directions these are maximal, we consider cross-sections of / = 0.5 
surface by equatorial plane {z — 0), by planes containing z axis and forming angle 
7r/4 or 7r/8 with xz plane, and finally by 2 planes parallel to xy plane, z = R/2 and 
z = R\/3/2. We mark results for these planes with index 0, 4, 8, 2 and 3 respectively. 
For each plane we consider a slice, subset of points on the surface / = 0.5 that are 
at most R/10 away from the crossing plane. For points on each slice we compute fit 
parameters ro and ri similarly to what we did for 2D flames. For instance, for plane 



2, parallel to xy plane, we flt cylindrical coordinates (pj = \J xf + yf, 4>i) of points 
in the corresponding nearby slice with pi — ro2 + ri2cos40^. 

The ratios ri/rg for the slices described are listed in Tab. 4.3 for models B and 
sKPP for several densities, for several moments of time for each density; these ratios 
give an idea about flame surface geometry. Moments of time shown were chosen 
as in the previous section. We remind that, according to 2D comparisons, open 
boundaries start influencing results for flame surface, for 256"^ runs, when its radius 



sKPP. 
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reaches approximately 120 km (about 60 cells) for central ignition, and approximately 
300 km for corner ignition (octant simulations). 
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Table 4.3: Aspliericity parameters, multiplied by 1000, at different times (indicated 
by radius R the flame has reached), for models B and sKPP (abbreviated to K in 
some places) in 3D. Flame velocity is D — 80kms~-'^, width l^i=3.2 cells for Model 
B, 4 cells for sKPP. Default initial flame radius is Rq — 30 km, box size 256^ (side 
528 km), central ignition by default. If other than default values were used, they are 
indicated next to the model name; octant simulations are indicated by letter "o" . 

As in 2D simulations, when expansion is small, 1/A < 0.4, both models behave 
well, asphericity parameters stay within 1% even when flame surface is close to box 
boundaries. According to theoretical estimate in Markstein approximation, quoted in 
the previous section, first unstable modes start growing at fiame radius exceeding 200 
cells at such expansions (for Model B), and growth rate is small at small expansion. 
Already at 1/A = 0.734, corresponding to fuel density 10^ g cm""^ in SNIa problem 
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(with half carbon, half oxygen fuel) asphericity parameters grow almost linearly with 
time for both models, although even extremal Ar/R remain within 5% for model B 
(as they do for this model for 1/A — 1.44). One can observe that rn is positive and 
grows for all times shown for model B, in contrast to the behavior for 2D simulations. 
It is noteworthy that for Model B results for octant and full-cube simulations are quite 
close to each other, with somewhat worse agreement for cross-section 3 (explicable, 
as flame surface forms sharp angle with crossing plane in this section, thus small 
3D perturbations would look larger in the section). This implies that reflecting 
boundaries do not alter signiflcantly LD instability of adjacent flame surface, thus 
octant simulations may be reliably used for modeling SNe la with central ignition 
(as far as "flame" propagation anisotropy and LD instability are concerned). 



CHAPTER 5 
MARKSTEIN EFFECT 



In numerical simulations of flame propagation in 2D and 3D we observed that the in- 
stantaneous fiame speed was somewhat different from the one measured for ID flame 
in a tube with the same diffusivity and burning rate in governing equation (1.3). It 
was smaller at the beginning of the simulation (at least for small expansion); Mark- 
stein effect is a possible physical explanation of this phenomenon. Flame front is 
stretched by geometrical effects when it is not planar; density and velocity distribu- 
tion across such a curved flame differ from those in planar steady flame, as a result so 
does instantaneous normal propagation speed. On general grounds one might expect 
corrections to flame propagation speed for small curvatures as expansion 



for the speed Dr^ of a flame with radius of curvature tq in inverse powers of (large) tq. 
In 3D this expansion could be more involved, with two principal radii of curvature of 
the flame surface tqi, ro2 appearing on the right hand side. Such systematic deviation 
of the speed of the curved flame from planar one was observed in experiments (see 
Markstein (1964)); the leading term, Ma/vQ is enough to account for curvature 
corrections in those. In 3D setting l/rg = I/tqi + 1/to2 must be used in the term 
linear in curvature (the only linear combination invariant under rotations). 

Negative sign at Ma/ R in (5.1) was chosen for convenience: with such a conven- 
tion Ma > for our model flames; some authors use notation with all pluses in (5.1). 
Notation Ma2 (as well as Ma^, and so on for coefficients of next terms in expansion 
(5.1)) is ours, such corrections of higher order in curvature are not studied in the 
literature. 

In this chapter we study this curvature effect, first in quasi-steady state approach 
we describe, and then numerically for the flame models presented in the previous 
chapter. We compare the results obtained using these 2 methods, and conclude 
that for model B physical Markstein effect dominates numerical effects, as well as 
corrections to effective flame propagation speed due to hydrodynamic instability 
development at expansions of 1/A = 1.44 and below. 



Here we present a general method for obtaining Markstein lengths of diffusion- 
reaction flames. No assumptions of small reaction zone thickness are made (unlike 
existing studies; for the model flames we consider, and are interested to apply the 
technique for, "preheating" and "reaction" zones have the same spatial scale). For 





5.1 Quasi-steady state technique 
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a spherical flame corresponding quasi-steady velocity (slowly changing with time as 
the flame radius changes) is found as an eigenvalue of a set of differential equations 
for quasi-steady distribution of matter velocity and quantities governing reaction 
rates (species concentrations and the temperature for physical combustion systems, 
reaction progress variable / for the model flames we considered before). Results for 
model flames are presented in the next sections. 

For clarity we present the method for a system with reaction rate dependent 
on concentration of one species / and on temperature T (a typical approximation 
in studying premixed flames in chemical and astrophysical simulations). For con- 
venience we suppose reaction rate $ = -R$o represented as a function of enthalpy 
density H instead of the temperature. Neglecting viscosity the system is described 
by Euler equation for matter velocity u, 

^ + (uV)u = -^+g, (5.2) 

mass continuity equation for evolution of the density p distribution 

^ + V(pu) = (5.3) 

and the following equations describing species diffusion and thermal conduction: 

1^ + uV/ = V{KVf) + i?$o(/, H) (5.4) 
dH 

— + uVH = V{n,VH) + gi?*o(/, H). (5.5) 

Here q represents speciflc heat release of the reaction. K = KKq^J, H) is diffusivity 
of the species described by variable f,K — K,KQ{f, H) is heat diffusivity; form- factors 
Kq, kq and $o ^-^e dimensionless. We assume / normalized as for the model flames 
of the previous chapters, namely / = corresponding to unburned fuel, and / = 1 
to reaction product (say, a state with dcflcicnt rcactant completely depleted). No 
external volume force pg is supposed to be present below. We consider isobaric 
burning regime (quasi-steady deflagration with matter velocities much smaller than 
the sound speed), as we did in Chap. 2. 

For model flames with reaction rate dependent on reaction variable / only the 
system is simpliflcd, in the same way as it is simplifled in the case of unit Lewis 
number Le (see below); essentially, Eq. (5.5) decouples from the rest of the system. 
With several species taken into account in reaction rates the needed modifications 
are also straightforward: an extra equation of type (5.4) is added to the system 
for each additional species; to the final eigenvalue problem (see below, (5.7-5.9) 2 
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ordinary first-order equations are added for eacli sucli species, no new complications 
are introduced. 

The system of combustion equations has a solution describing a spherical flame 
propagating outwards. It approximately describes deflagration in a system consisting 
of pure fuel at rest at t = and ignited at one point, after certain quasi-steady distri- 
bution of matter velocity and thermodynamic quantities is established, on timescales 
comparable to sound crossing time of the system. Such centrally symmetric distribu- 
tions (/(r), etc.) propagate with the flame, changing slowly due to geometry (flame 
radius) changing — in contrast to a planar flame there is no exact translational sym- 
metry (in r) for such a system, thus no steady solutions in the strict sense. We refer 
to system behavior like this as quasi-steady spherical flame propagation — when all 
quantities depend only on the distance from the flame center and time, and their 
changing with time in the frame coexpanding with the flame, that is expanding with 
so defined instantaneous flame velocity D{t), is much slower than in the rest- frame. 
For any quantity (/, p, H, radial matter velocity Ur, etc.) in coexpanding coordinates 
(r, t) {rjjt), 7] = r — J D{t)dt partial time derivatives are small at fixed location 
T) with respect to expanding fiame surface. For corresponding time derivatives in 
laboratory coordinates (r, t) 



dt 



dl 
dt 



(5.6) 



that is the changes due to profiles slowly reorganizing with changing geometry are 
negligible in comparison with changes due to fiame propagation. We rewrite the sys- 
tem (5.2-5.5) below in comoving quasisteady coordinate 7] thus getting rid of explicit 
time dependence (it will enter only implicitly through fiame radius rg dependence on 
time) based on approximate equality like (5.6) holding for p, i/, Ur- Such a substitu- 
tion for time derivatives yields the eigenproblem, which determines the quasi-steady 
fiame propagation speed D we seek for. 
In dimensionless variables 



fj = 


v^R/k 


f — 


r^R/k 


u — 


Ur/VRk 


d = 


D/VRk 


H = 


H/q 


A = 


k/k 
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after quasi-steady substitution of the form (5.6) equations (5.3-5.5) are rewritten as 

(u - d)^ = -pf-^^(ur^) (5.7) 
drj drj 

^n-d)^ = ^!^ + ^^ + Mf) (5.8) 

dfj drj drj r drj 



[u — d)—— = A 
drj 



d KQdH ^ anQ dH 



drj drj r drj 



+ Mf)- (5-9) 



In the above a + 1 — dimension of the problem: a = 1 for a cyhndrical flame, a — 2 
for a spherical flame (jS^) in 3D. Using corresponding equation of state and assumed 
isobaricity one expresses p as a function of H in Eq. (5.7), thus making the system 
closed. We used p = const/if as we did in the previous chapters. 
Boundary conditions for this problem we use are 

f = fo = ro^R/K : / = 1, df/dfj = 0, dH/dfj = 0, = 0; (5.10) 
r^+^: f^o,H^Ho/q = //(Pf,eb %el)/9- (5-11) 

It is through these boundary conditions that the flame radius tq enters the eigen- 
problem, thus yielding its eigenvalue d dependent on rg. In prescribing boundary 
condition at r = ro we assumed that the flame does not have inflnite tail into ash. 
Formulation of quasi-steady problem for systems with such an infinite tail encounters 
several problems, we will not touch them here. 

As in ID case the system simplifies when k, = K {i.e. Lewis number equals 1; we 
see, below, that these equal transfer coefficients may still depend arbitrarily on / and 
H for the simplification to hold; Le ~ 1 in terrestrial fiames, in nearly ideal gases). In 
this case we just reproduce a known result (Zeldovich & Prank-Kamenetskii (1938)) 
that distributions of / and H is similar — for any flame geometry (not necessarily 
spherical): Eq. (5.5) in this case coincides with (5.4) after rescaling H by q, and 
boundary conditions (for / and H going to constants behind and in front of the 
flame) then yield algebraic relation H — Hq + qf; one is left with a reduced system 
of equations (5.3-5.4). The same simpliflcation occurs in the case of model flames 
studied for Flame Capturing applications, when H — Hq + qf is postulated by the 
technique. 

We present results for D^tq) in the next section, for model flames studied in the 
Chap. 2. 
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5.2 Quasi-steady results for models used in Flame 

Capturing 

Below we stick to a reduced form of the eigenvalue problem, with H — Hq + qf. 
pH = const at constant pressure is assumed. Although the system (5.7-5.8) is not 
tranlationally invariant in f (unlike the case for a planar front) due to f explicitly 
appearing, f does not change significantly within the flame width at large radii^. For 
numerical solution we found it convenient to integrate over / instead of f, similarly 
to what we did in Sec. 2.1.2 in ID. This is especially successful for flame models with 
flnite total proflle widths. For such models we integrate the following system: 

I = aKoim-'/^ - ^^^^ + . (l + /) (l + (5.12) 

I ^ (5.13) 
df ^ n 

I _,.-V«^MZ)_l_. (,14) 



Here 



uq is a constant such that / — >-Oat/— >0(r— > +oo). Eq. (5.14) is (5.7) rewritten 
in terms of /, (5.13) is a rewritten deflnition of n = —KqIJ) df /dfj, and (5.12) is 
Eq. (5.8) in these new variables. 

Boundary conditions at / = 1 (left boundary of the flame, f = fo) are 

/ = 1 : ^ = r^,n = 0. (5.15) 

At / = (for finite flames this is the right boundary of the flame, located at finite 
radius fQ + Wc, for flames with an inflnite tail into fuel f = oo at / = 0) the boundary 
conditions are 

/ = 0: 7 = 0, n = 0. (5.16) 

This system was integrated from / = 1 to / = 0. Values of Iq = I{f = 1) and 
d were estimated based on ID results (initial seeding), and then solved for exactly 
using Newton- Raphson algorithm to satisfy the 2 boundary conditions at / = 0. 



1. We have tried estimating Ma in the approximation of large fo, treating terms ~ 1/^ with 
positive powers in Eq. 5.12 as perturbations, using ID flame profiles /(r) in estimating such terms. 
Wc only present exact treatment of the problem here for brevity, as approximate consideration is 
not significantly simpler. Very rough fully analytical estimates may be obtained assuming, say, 
linear flame profile in perturbation terms; this leads to about 20% error in Ma for finite-width 
flames. 
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System (5.12-5.14) is singular at / = 1 as n(l) = 0. Because of this integration 
was actually started at I — fe for certain small positive /e (0.01 of bulk integration 
step); asymptotic expansions for ^, n, I were used as initial values at / = 1 — /e. For 
model B (Eq. 3.4) these are 



C = ko + — 



pra-1 



00 



l-ra)(l+A)' 



d 



1 - c(A) 



1 



1+* 



(5.17) 



For models A (Eq. 3.1), original model of Khokhlov (1995) (same expression for 
*^o(/)i ^ = in the expression for diffusivity K) and sKPP (3.3) initial asymptotic 
values used are 



n — 
e = 
/ 



\ 2 J 



(5.18) 



h - 2/, 



= V2$o(l) 



1/2 ^ 

^o(l + A) L 



\ i/Q'^O 2 / J 

for A, Khokhlov (1995) 



2ea(l-e/) for sKPP 



1 



+ — 
0' ^0, 



Resulting curves for (i(ro) are shown in Fig. 5.1 for Model A, and in Fig. 5.1 
for Model B, for 4 expansions each. d{rQ) / d{rQ = oo) is plotted against wi/tq, the 
slope of the curves at wi/tq = thus gives Markstein number M. Deviations from 
linear Markstein law can be seen at larger curvatures. Values of Markstein numbers 
are shown in Tab. 5.1 in the next section, compared to these numbers estimated in 
direct numerical simulations. 
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Figure 5.1: Dependence of flame speed d on curvature, presented in units of inverse 
flame width w\ for Model A. = d{rQ = oo). 
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Figure 5.2: Dependence d{wi/rQ) for Model B. 
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5.3 Markstein numbers from direct numerical 2D 

simulations 

In this section we estimate how observed numerical flame speed depends on flame 
curvature, whether one can consistently get Markstein length parameter Ma that 
would describe this dependence via (5.1) for a range of rg's as the flame grows in 
radius. If this is the case, one may further try to correct for this curvature effect 
on flame speed in multidimensional simulations (by adjusting diffusivity and burning 
rate parameters based on local flame curvature, so that to correct for the discrepancy 
in flame speed between plane and curved flames). Comparison between analytical and 
numerical results clarifies whether real physical Markstein effect dominates unrelated 
effects (say, growing asphericity of flame surface studied in the previous section leads 
to apparent growth of D estimated assuming spherical flame shape). 

The dependence on rg we observe deviates from hnear form ((5.1) with Ma2/Fp' 
and further terms in expansion omitted). Contrary to chemical flames, with widths 
much smaller than radii of curvature studied, our artificial "fiame" widths are much 
wider: Rq = 30 km we typically use as initial flame radius is just about 15 grid 
cells, and "flame" width Wi is about 4 cells. This difference makes such deviations 
expectable, after all Markstein "law" is just an expansion over small parameter, which 
is a ratio of flame width and flame radius of curvature (it is this ratio that determines 
how strong a stretch on the flame scale is). To make comparisons more consistent 
we use fits of observed fiame speed with (5.1) with 3 terms left in parentheses (that 
is, (5.1) with ". . ." omitted), with Ma, Ma2 considered free parameters. We use the 
same fits to get Ma, Ma2 for found analytically. Markstein numbers 

M = Ma/Wi 

are reported below (as we saw in Sec. 5.1 analytically, Markstein length Ma scales 
proportionally to flame width if one changes the latter by varying scale parameters 
K.R of the model while keeping diffusivity and burning rate form-factors KQ[f), 
$o(/) invariant), see Tab. 5.1. 

For each model and expansion parameter we present Markstein number Mgt found 
analytically, together with that found numerically, for several timesteps and initial 
setups used. Analytical results shown were obtained by fltting Dj-q found through 
quasi-steady technique presented in Sec. 5.1-5.2 with 3-term expansion of the form 
(5.1) for a set of approximately equidistant fiame curvatures l/rg, ranging from 1/400 

to 1/50 (dimensionless, rg is scaled with the same factor \J K / R as dimensionless 
coordinate x and width wi in Sec. 5.1). As a reminder for the sense of scale, fiame 
width is ~ 2 for Model B, and ~ 8 for sKPP for small expansion, see Fig. 3.2 
and its caption. This way, when fits are used with at least quadratic in curvature 
terms left in expansion, exact value found for Markstein number M is not sensitive to 
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precise set of radii used for the fit. Say, if more curvatures are added to cover a range 
of to 1/15, the M found for Model B at 1/A = 0.167 changes from 0.2145 to 0.2117 
(1.3%). Besides, when more terms are used in fit (higher powers in curvature left in 
expansion (5.1)) results for M do not change drastically as well. For the example 
above leaving terms up to cubic leads to M = 0.2148 (for curvature range to 1/15; 
even smaller difference for reduced range l/rg G [1/400; 1/50], as for M's presented 
in Tab. 5.1); fit with terms up to quartic in curvature yields M = 0.2150. 

Numerical Markstein numbers were estimated from 2D runs with central and 
corner ignition (same setup as in the previous sections). For this, flame speeds and 
radii were recorded for each timestep; then D,.q was fitted with expansion of the form 
(5.1) with three terms (up to ) for timesteps from 2000 to the step when flame 
radius reached the value of R recorded in the column header in Tab. 5.1. Flame speed 
was estimated based on integral burning rate, D — — ^ ^ a ^"1^^^^ ; A here denotes 
flame surface, estimated based on burned volume assuming spherical flame shape. 



Model 






1/A 


M, i?2 


M, i?3 


M, i?4 


M, i?5 


Mst 


B 






0.167 


0.233 


0.232 


0.232 


0.232 


0.2145 


sKPPq, 


AT = 512 




0.167 


0.0360 


0.1506 


0.1106 


0.1107 


0.0733 


Bq 






0.398 


0.249 


0.241 


0.241 


0.242 


0.2163 


sKPPq 






0.398 


0.169 


0.159 


0.159 


0.158 


0.0680 


Bq, N = 


= 512 




0.734 


0.325 


0.309 


0.313 


0.326 


0.3246 


sKPP 






0.734 


0.233 


0.280 


0.409 


0.587 


0.0588 


sKPPq, 


TV = 512 




0.734 


0.251 


0.294 


0.416 


0.584 


0.0588 


B 






1.44 


0.384 


0.309 


0.305 


0.325 


0.4409 


B,i?o = 


15 




1.44 


0.202 


0.255 


0.258 


0.286 


0.4409 


sKPP 






1.44 


0.705 


2.113 


3.633 


6.843 


0.0517 


B,Ro = 


90, N = 


2ll 


2.33 




0.629 


1.094 


2.28 


0.3432 


B,Ro = 


90 




2.33 


-1.69 


0.760 


1.134 


2.48 


0.3432 


B,i?o = 


90, Wi = 


= 4.5 


2.33 


-1.13 


0.512 


0.364 


0.430 


0.3432 


B,i?o = 


90, Wi = 


-- 5 


2.33 


-1.64 


0.517 


0.325 


0.339 


0.3432 



Table 5.1: Markstein numbers computed numerically for models B and sKPP at 
timesteps, computed when flame radius reached R2 — Rq + 90 (km), R^ — Rq + 
270, R4 = Rq + 420 and R^ = Rq + 600 (the radius is shown as an index in the 
table header). Last column shows Markstein number found via quasi-steady state 
technique. 

As follows from the Table, (5.1) consistently describes dependence of flame speed 
on its radius at small expansions 1/A, one gets almost the same Markstein numbers 
while fitting Dji over different ensembles of -R's, meaning that (5.1) is a good ap- 
proximation, omitted terms are not significant for the curvatures studied. As fiame 
remains circular at these expansions for small enough radii, one would not expect 
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explicit dependence of flame speed on time, all time dependence is contained in Dj^(^£^ 
dependence, through Markstein effect and radius growth with time. Deviations from 
this are observed at larger expansions, again expectedly, as corrections to flame speed 
due to asphericity development (both large-scale, and LD-hke) depend exphcitly on 
time provided for the asphericity to develop. For sKPP these deviations are pro- 
nounced at 1/A = 0.734 and above, in accord with our direct observations of flame 
distortion development in the previous section. Model B shows similar violation of 
(5.1) in simulations at 1/A = 2.33 (pfuel = 10^ g cm~"^ in SN la problem). This 
again is in agreement with worse model behavior (faster LD-instability development) 
at such expansion. Notice that precise M found numerically differs from that found 
through steady-state technique increasingly as expansion increases. Same increasing 
deviation was observed above for D in ID simulations, thus it can be expected that 
similar discrepancy should be found for curvature corrections to D. Increasing flame 
width damps LD instability development, as well as decreases discrepancy (due to 
discretization) between numerically found Markstein number, and steady-state one 
(in continuous calculation). This can be seen in the Table for Model B entries at 
1/A = 2.33 and different flame width. Notice also, that for R = R2 aX 1/A = 2.33 
one gets negative Markstein numbers. Burning is not quasi-steady yet at that time 
(timestep 1300) for larger expansions. Studying flame speed as a function of time 
(what we effectively do through studying Markstein number evolution) provides an- 
other way to characterize flame asphericity development. 

For smaller expansions agreement with steady-state results for Model B is en- 
couraging. Similar agreement, with further deviation in the same direction due to 
discretization effects, is observed for 3D simulations; say, we get M — 0.280 for Model 
B 256^ cubic run at 1/A = 0.398 R — R2 — 120 km (curvature of a spherical flame 
in 3D being 2/R). 



CHAPTER 6 
CONCLUSIONS 

We analyzed a reaction-diffusion system (1.3) with hydrodynamically determined 
advcction velocity u. This model describes premixed physical flames with one dom- 
inating reaction, the rate of which can be effectively determined by one reaction 
progress variable /, which is the case, in part, for systems with unit Lewis number 
(ideal gases in part); this system is used as a tool for tracking unresolved flames 
in deflagration simulations (Flame Capturing technique, FC, Khokhlov (1995)), in 
part in astrophysical simulations of nuclear deflagration in a White Dwarf during an 
SN la explosion. 

One part of the study was performed in continuous steady-state approximation 
(steadily propagating flame solutions) for planar flames (spatially one- dimensional 
problem; Chap. 1 and 2) with the purpose of flnding the flame propagation speed and 
width, as well as for spherical flames in 2D and 3D, with the goal to quantify Mark- 
stein effect, how the flame front speed depends on the its curvature (Chap. 5). For 
this part of the study we assumed heat release proportional to / increase, 5Q = qdf 
{q meaning total heat release in the reaction, per unit mass), negligible pressure jump 
across the flame (thus, in part, thermal enthalpy density increasing linearly with /, 
dH — SQ), pH — const in the process of such isobaric burning. The problem of flnd- 
ing steady flame profiles f{x) {f{r) for spherical fiames) was analytically reduced 
to finding eigenfunctions of boundary- value problems, (2.10) with p{0) = p{l) = 
for planar flames, (5.12-5.16) for spherical flames. These eigenfunctions yield flame 
proflles in dimensionless spatial coordinate: all the quantities were made dimension- 
less by scaling with appropriate powers of characteristic values of reaction rate R 
and /-diffusivity K, Eq. (2.7). Dimensionless flame speed d, Eq. (2.9) is given by an 
eigenvalue of the corresponding boundary problem in ID or nD; for the case of spher- 
ical flames flame radius rg enters explicitly the boundary conditions, thus making 
flame speed and proflle dependent on flame curvature. 

Another part of the study was performed numerically, using full hydro dynamical 
codes ALLA (Khokhlov (1998)) and FLASH (FryxeU et al. (2000)). This part had two 
major goals. The first one was to check how based on steady analysis calibration of 
different flame models studied was affected by discretization effects in direct simula- 
tions. For flame widths used/proposed to be used in FC (3.2-4 cells between / = 0.1 
and / = 0.9) flame velocities observed in direct simulations are somewhat smaller 
than their prescribed values based on steady-state analysis, the difference increases 
with parameter 1/A characterizing matter expansion, Eq. (2.4); this difference re- 
mains within 8% for the models studied at 1/A < 1.44 (interval of most interest to 
SN la problem, corresponds to fuel (with composition 0.5 -"^^C -|- 0.5 -"^^O) density of 
3 X 10^ gcm~"^ and above). Flame widths observed in these simulations are larger 
than their prescribed values; Cf. Tab. 2.2 and 3.1 for comparison of continuous and 
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discrctizcd speeds and widths for some models studied. The discrepancy in width 
was clearly related to numerical diffusion tails (Fig. 2.9) by increasing number of 
computation grid cells within the flame width, Tab. 2.3; for such wider flames front 
speed also tends to its value prescribed through continuous steady-state analysis, 
discrepancy reduced to about 1% for 16 cells within flame width Wi. 

For thin flames used for FC it is desirable to correct for this discretization dis- 
crepancy; this is accomplished by deflning dimensionless flame speed d and width w 
by direct simulation in ID, using flame parameters ensuring physical flame width W 
the same as the one to be used for actual simulations (with gravity, arbitrary flame 
geometry, etc.), see Sec. 3.1.2 for more detail. 

For values of d and w found by any of the methods (for any flame model, specified 
by dimensionless form-factors of reaction rate and diffusivity $o and Kq, (2.7)) the 
procedure to obtain proper normalization factors R and K so that Eq. (1.3) coupled 
to the rest of hydrodynamic equations yield a "flame" having prescribed width W and 
speed D is given by Eq. (2.37) as described in Sec. 2.5.2. For Model B, in particular, 
d and w are given by flts (A. 6), thus yielding efficient procedure for prescribing K 
and $ based on local A. 

Any model with unique propagation speed may be calibrated this way to yield 
a flame with required properties, in idealized conditions at least, in ID with sim- 
ple hydrodynamics. Some models are however better than others for use in flame 
capturing. Certain features of $o a'^d Kq of the model used produce flames with 
undesirable features. Some of these features can be identified in steady state study: 
for fiamc models with Kq = 1, the original one (Khokhlov (1995), Eq. (2.11)) and 
KPP (2.12) $0 goes to zero either at / = or 1, thus producing flames with inflnite 
tails. For artiflcially thickened "flame" used in FC to deviate least in behavior from 
much thinner physical flame it is imperative to have the flame completely localized 
in the narrowest possible region. Central region of the "flame", characterized by 
large gradients of / must be adequately resolved in simulations, thus have physical 
width of about 3 or more grid spacings; therefore flames without long tails, with 
proflle close to linear are to be preferred. By small perturbation of the burning rate 
model KPP may be transformed into a model with flnite tails (and unique eigenvalue 
for propagation speed), sKPP, Eq. (3.3). However unless the perturbation (given by 
shift parameters ea, e/) is large enough (leading to significant numerical noises) model 
sKPP retains long tails in the profile (Fig. 3.4), thus having a disadvantage compared 
to other models proposed. 

Using non-constant diffusivity, Kq = f" allows one to make a model with step- 
function burning rate have finite tails; choice of parameters {fQ\r) = (0.2; 0.75) 
makes fiame profiles insensitive to fiame expansion, another desirable feature to make 
"flame" response to hydrodynamics similar at different densities (in SN la problem), 
consistent throughout the simulation as the flame propagates into less dense regions 
further from the center of the WD. However this model. Model A, as well as other 
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models with large discontinuities in burning rate produces significant numerical noises 
in simulations, about 10% fluctuations in instant flame speed, Fig. 3.5. Numerical 
noises, as well as flame propagation anisotropy, small-scale hydro dynamical insta- 
bility of the flame front in 2D and 3D are only observed in actual simulations of 

discrctizcd on computation grid set of hydrodynamic equations, not in continuous 
steady-state study. Studying these intrinsically non-steady features of the flames was 
another major goal of numerical part of our study, presented mostly in Chap. 3 and 
4. 

Numerical study showed that different models differ one from another in amount 

of noise they produce (which is readily observed in ID), in flame surface distortion 
in 2D and 3D. ID noises and flame speed anisotropy with respect to the grid are 
numerical artifacts and must be minimized by choosing appropriate flame model 
for use in FC. Small scale instability we observed on spherical, cylindrical and per- 
turbed planar flames is of physical nature, characteristic for any flame propagation 
when fuel and ash densities differ. This instability is of the type studied by Landau 
and Darrieus in simplified setting, as was demonstrated by quantitative agreement 
between characteristic instability lengthscales in our direct simulations, and analytic 
estimates in Markstein approximation. Sec. 4.4. It is physically impossible to invent 
a flame model immune to this instability. However some models, hke sKPP, become 
fast perturbed on very small scales, of order 10 cells at 1/A > 0.734, whereas another 
model proposed. Model B, Eq. (3.4-3.7), only starts showing signs of this instability 
development (at 1/A > 1.44) when getting close to boundaries in 1024 x 1024 box, 
in cylindrical flame simulations; flame radius increasing by a factor of 20 by that 
time. Fig. 4.1. Such behavior is adequate for SN la simulations. Numerical effects, 
propagation anisotropy and noises, arc also insignificant for Model B, whereas sKPP 
shows significant anisotropy at 1/A > 0.734, which strongly distorts the pattern of 
LD-type instability. 

Flames of Model B are not perturbed significantly by reflecting boundaries, in 
contrast with sKPP, for which octant simulations are thus even more question- 
able. Markstein numbers computed for Model B using quasi-steady state technique. 
Sec. 4.4.1 and estimated based on radius dependence of cylindrical fiame speed -D(ro) 
in simulations demonstrate close agreement at small expansions, and are close enough 
at 1/A = 1.44 to suggest that physical Markstein effect dominates over contribution 
to D{rQ) dependence due to propagation anisotropy and LD- instability at densities of 
> 3 X lO^g cm~^ in WD problem. Based on all these features observed we recommend 
using Model B for Flame Capturing. 

A few directions to develop that are close to our study are the following. 
• Larger 3D simulations need to be performed to clearly see LD-instability on 
spherical fiames for Model B. The results might suggest slight changes in the model, 
better behavior in 3D may be possible with somewhat larger c(A) term in burning 
rate at large expansions. 
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• To use the model when Hp is significantly nonconstant (which is the case in 
SN la problem at smaller densities, ~ lO^g cm~"^, near quenching of nuclear burning) 
steady-state model calibration has to be modified; this is automatically corrected for 
if the model is calibrated numerically, as is the case for Model B calibration (A. 6). 

• Also at smaller densities in SN la problem assumptions of the fiamc being 
isobaric becomes increasingly violated. This is due to decreasing sound speed, and 
increasing speed of the fiame with respect to the grid, equal to + 1/A) for planar 
fiames. As follows from analysis in Sec. 2.1.1 fiame speed depends on exact density 
distribution within the fiame in part, which is affected by pressure jump across the 
fiame and should be corrected for if this jump is significant. The jump depends on 
geometry of the fiow: in our simulations this non-isobaricity led to systematic in- 
crease of fiame speed in 2D and 3D as compared to the speed in ID. The deviation 
reached 1% at 1/A = 2.33, and ~ 20% at 1/A = 8.57 (at that density, 8 x lO^g cm'^, 
sound speed is just a factor of 2 larger than D{1 -\- 1/A), for D = 80 km s"-*^ used.) 
These problems of equation of state and large pressure jump should not be relevant 
for most terrestrial fiames. 

• Quasi-steady state technique presented in Sec. 5.1 in general setting, with 
arbitrary Lewis number and possibility of having several species involved in reac- 
tions with comparable rates is immediately applicable for quasi-steady estimates of 
Markstein widths of various terrestrial and astrophysical problems. The same tech- 
nique, simplified to planar fiames, may be used to find steady propagation velocities 
of fiames more general than the ones described by Eq. (1.3). This may be needed 
even for some FC models that involve more than one reaction progress variable, with 
intent to more accurately prescribe heat release distribution when several character- 
istic reactions with different enough characteristic temporal and spacial scales are to 
be taken into account. 



APPENDIX A 
SUMMARY OF IMPLEMENTATION OF FC 
TECHNIQUE WITH FLAME MODEL B 

Here we summarize a step-by-step instructions for using the Flame Capturing Tech- 
nique based on Model B we recommend based on analysis in the thesis. 

Assuming one has a hydrodynamic solver for modeling deflagrations, the steps to 
use the technique for tracking the thickened flame region are the following: 

1) Add one new scalar quantity / to physical variables evolved with the code. 

2) Use 5Q/dt — q df / dt as a heat release term in original hydrodynamic equations. 
q is total specific heat release of nuclear burning, it depends on local pressure and 
fuel composition. 

3) Evolve / via 

^ + u-V/ = V(i^V/) + *(/), (A.l) 
where u stands for local gas velocity, source term and diffusivity are given by 

$(/) = ^r^(i-/r'^(/-c(A)), (A.2) 

K{f)^Krf{l-fY-. (A.3) 

respectively, 

= Sa = 0.8, = 1.2, ra = 0.8. (A.4) 

Normalization factors K and R are functions of locally determined expansion param- 
eter A = Pash/ (Pfuel ~ Pash) (the latter depends on pressure and fuel composition in 
the cell), and the values for "flame" speed D and width Wi (between values / = 0.1 
and / = 0.9) one strives to obtain: 

d(A) wi(Ay d{A)/wi{Ay ^ ' 
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Dimensionless d{A), w(A), and term c(A) in Eq. (A. 2) are given by the following 
fits: 

= 0.005 

1/A e [0; 0.515] : { d= 0.328227 - 0.10051/A + 0.0244596/A2 
n = 2.075422 + 0.443918/A - 0.097483/A2 
= 1/A -0.51 

1/A e [0.515; 0.81] : { d ^ 0.497578 - 0.363476/A + O.IOSG/A^ 

71 = 1.553031 + 1.50381/A - 0.192923/A2 
= 0.3 

1/A e [0.81; 1.5] : { d= 0.172133 - 0.051673/A + 0.0073512/A2 

n = 2.475085 + 0.232926/A - 0.0323989/A2 
= 0.675 - 0.25/A 
1/A e [1.5; 1.9] : { d^ -0.0420695 + 0.129058/A - 0.0177843/A2 
n = 2.764056 - 0.0061156/A + 0.0025170/A2 
= 0.2 

1/A e [1.9; 8.6] : { d= 0.0139649 + 0.434752A - 0.507838A2 + 0.250103A3 
wi = 3.3891706 - 2.857323A + 5.045909A2 - 3.744233A3. 

(A.6) 

These fits yield errors in flame speed in ID not exceeding 0.8% for Wi = 3.2 cells, 
the value we recommend to use, and 1/A e [0;8.6], covering expansions of interest 
for SN la problem. 
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